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ABSTRACT 

Large surveys using the Sunyaev-Zel'dovich (SZ) effect to find clusters of galaxies are 
now starting to yield large numbers of systems out to high redshift, many of which 
are new discoveries. In order to provide theoretical interpretation for the release of 
the full SZ cluster samples over the next few years, we have exploited the large- 
r-j '. volume Millennium Gas cosmological iV-body hydrodynamics simulations to study 

O f the SZ cluster population at low and high redshift, for three models with varying gas 

physics. We confirm previous results using smaller samples that the intrinsic (spherical) 
Y500 — -M500 relation has very little scatter (ciog 10 y — 0.04), is insensitive to cluster gas 
physics and evolves to redshift one in accord with self-similar expectations. Our pre- 
heating and feedback models predict scaling relations that are in excellent agreement 
with the recent analysis from combined Planck and XMM-Newton data by the Planck 
Collaboration. This agreement is largely preserved when r^g and M500 are derived 
using the hydrostatic mass proxy, lx.500> albeit with significantly reduced scatter 
("log y — 0.02), a result that is due to the tight correlation between I500 and Ix.soo- 
Interestingly, this assumption also hides any bias in the relation due to dynamical 
activity. We also assess the importance of projection effects from large-scale structure 
along the line-of-sight, by extracting cluster I500 values from fifty simulated 5° x 5° 
sky maps. Once the (model-dependent) mean signal is subtracted from the maps we 
find that the integrated SZ signal is unbiased with respect to the underlying clusters, 
although the scatter in the (cylindrical) I500 — Afsoo relation increases in the pre- 
heating case, where a significant amount of energy was injected into the intergalactic 
medium at high redshift. Finally, we study the hot gas pressure profiles to investigate 
the origin of the SZ signal and find that the largest contribution comes from radii close 
to r5oo in all cases. The profiles themselves are well described by generalised Navarro, 
Frenk & White profiles but there is significant cluster-to-cluster scatter. In conclusion, 
our results support the notion that I500 is a robust mass proxy for use in cosmological 
analyses with clusters. 

Key words: hydrodynamics - methods: numerical - X-rays: galaxies: clusters - galax- 
ies: clusters: general 



1 INTRODUCTION 

The Sunyaev-Zel'dovich (SZ) effect (| Sunyaev fc Zel'dovichl 
1 19721 ) is a powerful method for discovering new clusters of 
galaxies. It arises generically due to the scattering of Cosmic 
Microwave Background (CMB) photons off free electrons, 
leading to a predictable spectral distortion in the CMB that 
is, in the non-relativistic limit, li nearly dependent o n the line 
integral of the electron pressure (|Birkinshawlll999l h In mod- 
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ern theories of structure formation, the dominant contribu- 
tion to the SZ signal comes from the intracluster medium 
(ICM), a diffuse plasma within clusters that is approxi- 
mately in hydrostatic e quilibrium within the dark-matter 
domi nated potential fsee lVoitll2005l : lAllen. Evrard fc Mantd 
l201ll for recent reviews). The key SZ observable is the Y 
parameter, defined as 



Y= / ydfi, (1) 
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where the integral is performed over the solid angle sub- 
tended by the cluster. The Compton-y parameter is deter- 
mined by the thermal structure of the ICM 



V 



n c T c dl, 



(2) 



where n c and T B are the density and temperature of the free 
electrons respectively and dl is the differential line element 
along the line-of-sight. Since Y can be expressed as a volume 
integral of the pressure (when the redshift and cosmologi- 
cal parameters are specified), it measures the total thermal 
energy of the gas, a property that ought to be strongly cor- 
related with the cluster's mass through the virial theorem. 
This means that Y ought to be relatively insensitive to the 
complex micro-physics taking place in the cluster core, un- 
like other global properties such as X-ray luminosity. 

Early observational studies confirmed the detection of 
an SZ signal towards known massive clusters of galaxies and 
used this to estimate the Hubble c onstant (e.g. I Jones et al.l 
1 19931 ; Birkinsh aw fc Hughes! fl994V Over the past decade, 
SZ observations of known bright X-ray bright clusters have 
become routine, allowing the inves tigation of cluster scal- 
ing relations t o be p erformed (e.g. iMcCarthv et al.l l2003b; 
Benson et"afl 12004 iMorandi. Ettori fc Moscardinil |2007t 



Bonamente et alT 2009; Huang et al.l 2010l ; Shimwell et al.l 
20111 ; lLancaster et alj 1201 J >. One potential shortcoming of 



this approach is that the samples are X-ray selected and 
therefore biased towards luminous, cool-core systems at low 
redshift. 

In the past few years, SZ science has entered the 
exciting new phase of blind surveys, where detections 
of ne w clusters have become possible (|Staniszewski et al.l 
2009). Indeed, SZ surveys are now yielding large num- 
bers of SZ-selected clusters, many of them new de- 
tections, especially from the South Pole Te l escope 



fSPT; IVanderlinde et ail |2010| ; lAndersson et alj l201ll ; 



IWilliamson et alj|201ll) , ~the Atacama Cosmology Telescope 



(ACT; iMarriage et all l201ll; [Sehgal et al.l l201lh and the 
Planck satellite l|Planck Collaboration! 2011a b e). Since the 
SZ effect is effectively independent of redshift, the SZ selec- 
tion function tends to favour higher redshift systems than 
the X-ray counterpart, assuming similar angular resolution. 
As a result, the new blind SZ surveys are starting to find 
new massive systems at z ~ 1 (Planck Coll aboration! [201 ldl : 
iFolev et al.ll201ll ; lMenanteau et alj201ll ). In the near future, 
we should expect to see these numbers increase substan- 
tially as the full survey results are published, nicely comple- 
menting X^ra^_s_urvejrs_sucli_asti Cluster Survey 
(MACS; lEbehng. Edge fc Henrvll200ll) and the XMM Clus - 
ter Survey fXCS; lRomer et al.ll200ll ; iMehrtens et alj|201ll ). 
Such complementarity will be further exploited with the 
next generation of X-ray surveys (e.g. with eROSITA) an d 
millimetric telescopes (e.g. CCAT; see lGolwala et ail [2009). 

One of the main goals of SZ surveys is to mea- 
sure cosmological parame t ers ( e.g. Barbos a et al. I 1 19961 ; 
ICarlstrom. Holder fc Reese! 12003 : IBattve fc Welled l2003h . 
Central to the cosmological application of SZ surveys 
is the scaling relation between the observables (Y and 
redshift, z) and cluster mass, M. Under the a ssump- 
tion that clusters form a self-similar population (|Kaiserl 
Il986l ) the SZ flux should scale as Y <x M 5/3 H(z) 2/3 , 
when measured within a radius enclosing a mean density 



that is a constant multiple of the critical density of the 
Universe. Early theoretical studies combined such sim- 
ple scaling relations wi th the Press-Schechter formalism 
dPress fc Schechterj[l97i) to predict the SZ evolution of the 
cluster popula t ion in a variety of cosmologi c al models (e.g. 
Cole fc Kaiseil ll98Sl; iBartlett fc Silkl Il994l; iBarbosa et all 



19961; lEke. Cole fc Frenkl Il996l; lAghanim et al.l 119971 ; 
Kav. Liddle fc Thomas! l200ll : IBattve fc Weller! |2"o03| ) . More 



recently, attention has turned to more detailed studies of 
how cluster gas physics impacts upon SZ scaling relations, 
both usi ng semi-analytic models ( e.g. IMcCarthv et al.l 
l2003al lbl: IShaw. Holder fc Bode! I2008T ) and full cosmolog- 



da Silva et al 



da Silva et al 

Bonaldi et all 12007 



2009 



ical A^-body/hydrodynamic simul a tions 
2000 [; IWhite. Hernauist fc Springell I2OO2I 
2004; iMotl et all 120051; iNagail 1200a ._ 
Hallman et al. 120071; lAghanim. da Silva fc Nunesl 
Battaglia et alj l201ll ). Simulations are now also be- 
ing used to invest igate the e f fects o f mergers on SZ 



scaling relations ilPoole et all 120071; IWik et alj 1200 



lYang. Bhattacharva fc Rickerl |2010| ; iKrause et alj l201ll 



A generic result from these studies is that the self-similar 
description appears to be approximately valid on cluster 
scales (M > 10 14 ft _1 MQ) but in detail, differences are 
seen between the models that are due to the effects of 
non-gravitational physics (cooling and heating processes), 
especially at low mass where the gas fraction is depleted. 

Two of the main shortcomings in previous simulation 
studies are the relatively small samples (that are some- 
times restricted to lower-mass clusters) and a limited range 
of (uncertain) cluster gas physics models, often not cal- 
ibrated to match X-ray data. Some studies may satisfy 
one of these criteria but usually not both. A new gen- 
eration of sim ulations are now s tarting to overcome both 
shortcomings. I&tanek et all (|20 10T ) recently pr esented results 
from two of the Millennium Gas simulations (|Hartlev et alj 
2008), large- volume runs based on the Millennium Simula- 
tion (ISpringel et alj|2005l ) with varying gas physics. These 
simulations are sufficiently large to enable the full range of 
cluster masses (10 14 - 10 15 fc _1 M©) to be studied and one 
of the runs, where the gas was pre-heated at high redshift, 
is able to match the mean X-ray lum inosity-temperature re- 
lation at z = ilHartlev et alJl20o"8T ). Although the work of 
IStanek et ail (|2010l ) was focused on the more general issue 
of multi-variate scaling relations, they presented results for 
the SZ Y — M relation measured within a radius correspond- 
ing to a mean internal density equal to 200 times the critical 
density, r 2 oo- 

The aim of this paper is to use these Millennium Gas 
simulations to focus in more detail on predictions of the SZ 
effect and, in parti cular, the Y — M rela tion for clusters. Our 
paper builds on the lStanek et al. (2010) work in three impor- 
tant ways. Firstly, we add a third model that includes a more 
realistic treatment of feedback, both from supernovae and 
active galactic nuclei. This model has already been shown 
to successfully m atch many of the X-ray properties of non - 
cool core clusters (|Short fc Thomas!l2009l ; IShort et alfeoiOT l. 
Secondly, we include in our analysis simulated maps of the 
full SZ effect along the line-of-sight, to assess the projection 
effects of large-scale structure. Finally, we attempt to pro- 
duce results for our Y — M scaling relations using methods 
that are more closely matched with observations. In partic- 
ular, we present our results for the smaller rsoo and investi- 
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gate the impact of assuming hydrostatic equilibrium and a 
mass proxy (Yx . iKravtsov. Vikhlinin fc Nagaill2006T ) on the 
Y — M relation. 

We organise the remainder of the paper as follows. In 
Section[2]we outline the simulation details and our methods 
used to define cluster properties. We also present some basic 
properties of the sample and SZ maps. Sections O [4] and [5] 
contain our main results: in Section[3]we present an analysis 
of the hot gas pressure profiles, before going on to study SZ 
scaling relations in Section [4] and the impact of hydrostatic 
bias in Section [5] Finally, in Section [6] we summarise our 
main conclusions and outline future work. 



2 SIMULATION DETAILS 

Our results are drawn from the Millennium Gas simulations, 
a set of large, cosmological hydrodynamics simulations of the 
ACDM cosmology (!T2 m = 0.25, Q. A = 0.75, fi b = 0.045, h = 
0.73, as = 0.9). In this section we summarise the details of 
these simulations and present our methods for constructing 
simulated cluster properties and SZ sky maps. 

2.1 Millennium Gas simulations 

The Millennium Gas simulations (hereafter MGS ; 
see lHartlev et al.l l200g|; Istanek. Rudd fc Evrardl 120091 ; 
IStanek et all l20ld ; IShort et al.l l20ld ; lYoung et ail l201ll 1 

were constructed to provide hydrodynamic versions of the 
Virgo Consortiu m's dark matt e r Mi llennium Simulation 
(hereafter MS; ISpringel et ai] I2005T ). The simulations 
were therefore started from the same realisation of the 
large-scale density field within the same comoving box-size, 
L = 500/i _1 Mpc and used the same set of cosmological 
parameters. The MGS were run with the p ublicly- availabl e 
GADGET2 A-body/hydrodynamics code (|Sprin gcl 2005). 
Due to the increased computational requirements from the 
inclusion of gas particles, the simulations were run with 
fewer (5 x 10 8 each of gas and dark matter) particles in 
total than the MS. The particle masses were therefore set to 
m gas = 3.1 x 10 9 /i _1 M and m dm = 1.4 x 10 10 fc -1 M s for 
the gas and dark matter respectively. Gravitational forces 
were softened at small separations using an equivalent 
Plummer softening length of e = 100 h~ kpc, held fixed in 
comoving co-ordinates. At low redshift (z < 3) the softening 
was then fixed to e = 25 /i _1 kpc in physical co-ordinates. 

Two versions of the MGS were run with the above prop- 
erties. Both runs started from identical initial conditions 
but differed in the way the gas was evolved. In the first 
run, the gas was modelled as an ideal non-radiative fluid. 
In addition to gravitational forces, the gas could undergo 
adiabatic changes in regions of non-zero pressure gradients, 
modelled using t he Smoothed Particle Hydrodynamics for- 
malism (SPH; see lSpringel fc Hernouistll2002l for the version 
of SPH used in GADGET2). Additionally, in regions where 
the flow was convergent the bulk kinetic energy of the gas 
is converted into internal energy using an artificial viscosity 
term; this is essential to capture shocks and thus generate 
quasi-hydrostatic atmospheres within virialise d dark mat- 
ter h aloes. In accord with previous studies (e.g. IShort et al.l 
120101 ) we refer to this simulation as the GO (Gravitation 
Only) model. 



It is well known that a non-radiative description of 
intracluster gas does not agree with the observed X-ray 
properties of clusters, especially at low masses, where an 
excess of core entropy is required to pro duce a ste eper 
X-ray luminosity-temperature relation (e.g. iVoitj [20051 ). A 
simple method capable of generating this excess entropy 
is to pre-heat the gas at high redshift before cluster col- 
lapse (|Kaisedll99ll ; lEvrard fc Henrvll99ll ). We implemented 
this method in a second simulation by raising the mini- 
mum entropy of the gas (by increasing its temperature) 
to -Kmm = 200 keV cm 2 at z = 4. The entropy level was 
chosen so as to matc h the mean z = X-ray luminosity- 
temperature relation (|Hartlev et al . 2008). We also included 
radiative cooling, an entropy sink. However, this made very 
little difference, as the cooling time of the pre-heated gas is 
very long compared to the Hubble time and therefore gas 
could no longer cool and form stars before the end of the 
simulation. We refer to this simulation using the label PC, 
for Pre-heating plus Cooling. 

We also consider a third model when analysing the 
SZ properties from individual clusters. Th is is the Feedback 
Only (hereafter FO) model developed bv lShort fc Thomas 

(2009) and then applied to MGS clusters by IShort et ~ 

(2010) , where full details of the method may be found. 
B riefly, it uses the semi-ana lytic galaxy formation model 
of |Pe Lucia fc Blaizotl (|2007T l. run on dark-matter-only res- 
imulations of MS clusters, to provide information on the 
effects of star formation and feedback on the intracluster 
gas. The model works as follows. Galaxy merger trees are 
first generated by applying the semi-analytic model to the 
dark matter distribution. Various properties of the galaxies 
(such as their position, stellar mass and black hole mass) 
are stored at each snapshot of the simulation. The clusters 
are then re-simulated with gas, assuming that the gas par- 
ticles have zero gravitational mass; this guarantees that the 
dark matter distributions (and therefore galaxy positions) 
are identical to those in the parent dark-matter-only sim- 
ulation. At each snapshot time, two important changes are 
made to the gas. Firstly the increase in stellar mass of each 
cluster galaxy is used to convert local intracluster gas into 
stars, a re quirement for gene rating sensible stellar and gas 
fractions (|Young et alj 1201 lh . This change in stellar mass 
is also used to heat the gas from supernova explosions. Sec- 
ondly, any increase in black hole mass is used to heat the gas 
on the basis that such accretion leads to an Active Galac- 
tic Nucleus (AGN) . The heating rate, known as AGN feed- 
back, is taken from lBower. McCarthy fc Benson! \200$) and 
is given by 

Ef ee d = min (esMBH^Edd, e r M B nc 2 ) , (3) 

where esmbh = 0.02 dictates the maximum heating rate (in 
units of the Eddington luminosity) and e r = 0.1 is the effi- 
ciency with which the accreted mass is converted into feed- 
back energy. This is particularly important because AGN 
are the dominant feedback mechanism on cluster scales. 

We ana l yse th e same sample of 337 clusters studied by 
IShort et all feoich . comprising all objects in the MS with 

1 In the usual way, we take entropy to mean the quantity K = 
kT/n^ 3 , where T is the gas temperature and n e the free electron 
density. 
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virial mass M V - 1T > 5 x 10 14 h x Mq and a random sample at 
lower mass (1.7 x 10 13 ft _1 M < M vir Ox 1O 14 /i _1 M ) 
chosen such that there were a fixed number of objects within 
each logarithmic mass bin. The FO model successfully gen- 
erates the required excess entropy of the low redshift popu- 
lation and provides a good match to the structural proper- 
ties of non-cool core clusters. The main shortcoming of this 
model is that it neglects the effects of radiative cooling and 
therefore cannot reproduce the m ost X-ray luminous cool 
core population |Short et al.l l201C0. This failure may not be 
as serious as it seems, however, since there is some evidence 
that the X-ray cool core population diminis hes with increas- 
ing r edshift, both from ob servations (e.g. jMaughan et aT] 
l201ll 'l and simulations (e.g. iKav et al.ll2007h . Furthermore, 
as we will demonstrate, the SZ Y parameter (which mea- 
sures the global thermal energy of the intracluster gas) is 
reasonably insensitive to changes in gas physics that pre- 
dominantly affect the cluster core. Issues relevant to this 
study where cooling could impact upon our results are the 
degree to which the ICM is hydrostatic and the effect of gas 
clumping on the X-ray quantity, Yx, used as a cluster mass 
proxy. We note that a first step towards including radiative 
cooling in the model has been made a nd shows promising 
results (|Short. Thomas fc Youn3l2012} ). Ultimately, a fully 
self-consistent scheme is desirable, where the same cooling 
and heating rates are used in both the semi-analytic model 
and hydrodynamic simulation. 



2.2 Cluster definitions and estimation of global 
properties 

Cluste rs are denned in exactly the same way as in lKav et al.l 
< |2007t ). Firstly, a friends-of-friends code is run on the dark 
matter particles for each snapshot. The dimensionless link- 
ing length (in units of the mean inter-particle separation) is 
set to b — 0.1, chosen to minimise the probability of linking 
two haloes together outside of their respective virial radii. 
The dark matter particle with the most negative gravita- 
tional potential energy is then identified for each group and 
this is taken to be the centre. 

In the next stage, a sphere is centred on each friends-of- 
friends group and its radius increased until the total mass 
(from dark matter, gas and stars, when present) satisfies 



M A 



4i j . , . 
— r A Ap CT (z) 



(4) 



where r& is the proper radius of the sphere, A is a specified 
density contrast, p C i{z) = {3Hq /8ttG)E(z) 2 is the critical 
density and E(z) 2 = Q m (l + z) 3 + J7a for a flat universe. 
We assume A = 500 for the main results in this study as 
this value is commonly adopted for observational studies 
(some of which we will compare to) because rsoo is suffi- 
ciently large to make many integrated properties insensi- 
tive to variations in core structure, while also being small 
enough to be within reach for detailed X-ray observations 
of many objects. We occasionally use the value of A appro- 
priate for the virial radius, r V ir, as defined by the spherical 
top-hat collapse model. This is a redshift-dependent quan- 
tity, A = A c (z), which we calculate using the fitting for- 
mula given bv lBrvan fc Normanl (|l998T ). Note that at z = 0, 
A c ~ 94 and r vir ~ 2r 5 oo- 

Once the cluster's mass and radius is defined, we cal- 



culate various properties of the hot gas, the most important 
being the SZ flux. The frequency independent part is given 
by 

*» = ^^/**z;dv; (5) 

where Da is the (cosmology-dependent) angular diameter 
distance to the cluster and the integral is performed over 
the entire cluster sphere. To simplify matters, we re-define 
the integrated SZ Y parameter 



Y 500 Dj 



Y500, 



(6) 



since this combination is directly proportional to the inte- 
grated thermal energy of the gas which is the physical prop- 
erty of interest. Note that the dimensions of Y500 are now 
that of area; we will therefore present values in h~ 2 Mpc 2 
units. The value of Y500 is estimated for each cluster using 



y 5 oo = 



fj, e mnm e c 2 



JVhot 

E 



Ti. 



(7) 



where the sum runs over all hot (T > 10 5 K) gas particles 
within rsoo, with mass m gas and temperature Tj. We adopt 
the value p c = 1.14 for the mean molecular weight per free 
electron, appropriate for a fully ionised plasma of hydrogen 
(with mass fraction X — 0.76) and helium (with mass frac- 
tion Y = 1 — X). We also assume equipartition of energy 
between the electrons and nuclei, thus T = T c . 

We estimate the X-ray temperature of the ICM us- 
ing t he spectroscopic-like temperature T s \ (Ma zzotta et al.l 
l2004h . appropriate for bremsstrahlung in hot (kT > 3keV) 
clusters 



EST PiT, 



1/4 



-3/4 ' 



(8) 



where p; is the density of particle i and in this case the 
sum runs over all hot gas particles with kTi > 0.5keV. We 
measure T s \ in the region outside the cluster core (a; corc < 
r/r$oo < 1, where x COIC = 0.1 for the GO and PC models, 
and 0.15 for the FO model ) to provide a closer match to 
observed X-ray temperature measurements (where a larger 
variation in core temperature is seen than in our simula- 
tions). 

A quantity related to Y500 is ^x,soo oc 
M eas Tx, estimated from X-ray data. Introduced by 
iKravtsov. Vikhlinin fc Nagail (|2006t ). it was shown to be a 
low-scatter proxy for cluster mass (due to scatter in X-ray 
temperature being negatively correlated with scatter in gas 
mass). We estimate this quantity as 



cr-rfc 



p c mum c c A 



1 T s i 



(9) 



where M gas ,5oo is the mass of hot gas within rsoo, although 
we occasionally present lx,soo in its native (/i _1 Mq keV) 
units, i.e. simply assuming ix.soo — Afgas.soo^Tsi- The main 

2 The GO/PC and FO data were processed independently and 
different choices for x COTC were made at those times. However, 
the effect of this difference on T s \ is small; we checked by re- 
calculating the GO/PC temperatures at z = using x corc = 0.15 
and found only a 2-3 per cent increase, on average. 
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Table 1. Number of clusters in our samples at redshifts, z = 0, 0.5 
and 1. Column 1 gives the model label and column 2 the redshift. 
Column 3 lists the total number of clusters in each sample with 
-^500 > 10 14 /i -1 Mq, while columns 4 and 5 sub-divide the sam- 
ple into the regular and disturbed populations respectively, using 
the s parameter defined in equation (1 1 D . 



Model 


Redshift 


N i 

1 Y clus 


iV c , us (s < 0.1) 


AT clu8 (s > 0.1) 


GO 


0.0 


1110 


986 


124 




0.5 


567 


457 


110 




1.0 


139 


103 


36 


PC 


0.0 


883 


799 


84 




0.5 


436 


355 


81 




1.0 


102 


78 


24 


FO 


0.0 


188 


154 


34 




0.5 


118 


122 


26 




1.0 


75 


51 


24 



difference between Y and Yx is that the former depends 
on the mass-weighted temperature while the latter depends 
on the X-ray tempera ture, which is more h eavily weighted 
by lower entropy gas l|Mazzotta et al.| [2004). Comparing Y 
with Yx therefore implicitly tests the dumpiness of the ICM 
since clumpy gas will be cooler and therefore lower the X- 
ray temperature rela tive to the mass-weighted temperature 
fe.g. lKav et~ai]|2008t ). As we show below, this effect is model 
dependent but is of minimal importance in the PC and FO 
simulations. 



2.3 Cluster sample 

Table [T] summarises the number of clusters in each of the 
runs at redshifts, z = 0,0.5 and 1. For our fiducial sample 
we have employed a lower mass cut of M500 > 10 14 /i _1 Mq, a 
useful limit for comparing with SZ cluster data. The GO and 
PC simulations have similar numbers, although the latter is 
slightly smaller due to the effect of pr e-heating on the gas 
fraction (|Stanek. Rudd fc Evrardll2009l ). Note the number of 
clusters at z = 1 is around an order of magnitude lower than 
at z = 0. There are significantly fewer FO clusters at any 
given redshift due to the fact that it is not a volume-limited 
sample. The drop in number at high redshift is not as severe 
in this case as the mean mass of the sample is higher and so 
a smaller fraction of clusters drops below the imposed mass 
limit. 

We also consider the effect of ongoing mergers by split- 
ting our sample into regular and disturbed sub-samples, us- 
ing a simple estimator known as the substructure statistic 
i|Thomas et al.lll998l ; iKav et al.ll2007l ). defined as 



JT.DO 



(10) 



where r<p is the position of the cluster centre (defined here 
to be the position of the dark matter particle with the most 
negative potential, <j>) and r cm is the centre-of-mass. We de- 
fine those clusters with s > 0.1 as disturbed systems and 
those with s ^ 0.1 as regular systems, although note that 
this terminology is strictly for convenience as all clusters 
are disturbed to some degree. In practice, this value delin- 
eates those that are clearly undergoing significant mergers, 



as discussed in IKav et al l |2007l ). The fraction of disturbed 
clusters increases with redshift in all models, from around 10 
per cent at z = to 25 per cent at z = 1, in the GO and PC 
models. Again, the different method for cluster selection in 
the FO model modifies the result but nevertheless the trend 
of increasing disturbed fraction with redshift is still seen. 

2.4 Cluster profiles 

We discuss hot gas pressure profiles in Section[3]as these are 
important for understanding the relative contribution to the 
SZ signal from different radii. The profiles are constructed 
by first identifying all hot gas particles within a radius rsoo 
of the cluster centre. This sphere is then sub-divided into 
spherical shells with fixed radial thickness in log 10 (a;), where 
x = r/rsoo. The pressure within the shell is then estimated 
using a mass-weighted average 



P{x) = 



k 



V 

V{x) fimu ^-^ 

i — l 



rmTi, 



(11) 



where the sum runs over all hot gas particles within the 
shell at radial position x, V is the volume of the shell and 
H — 0.59 is the mean molecular weight for an ionised plasma 
(assuming zero metallicity). 



2.5 Cluster maps 

We also compute the thermal SZ effect due to an individ- 
ual cluster by constructing Compton-t/ maps. This allows 
us to separate the cluster contribution (within a cylinder) 
from the total integrated signal along the line-of-sight. Each 
map is constructed by first identifying all hot gas particles 
within a cuboid of size 2r v i r x 2r- v i r x 6r v i r , centred on the 
cluster. The particles are then projected along the long axis 
of the cuboid and smoothed on to a 2D grid, creating the 
y distribution. We estimate y at the location of each pixel, 
R-p = (x,y), as 



w(\Rj - R p \,hi)Tj 

Api^ c mum B c 2 ^ V w(\Ki - R p |, hi) ' 

i — P 



(12) 



where A p i x is the area of a single pixel and w is the projected 
version of the SPH kernel used by GADGET2. The main sum 
runs over all hot gas particles with projected position R^, 
temperature T; and SPH smoothing length hi. The sum in 
the denominator runs over all pixels and normalises the ker- 
nel for each particle. 

Fig. [T] illustrates Compton-j/ maps for two massive clus- 
ters in our simulations at z = 0: a regular (s ~ 0.02) cluster 
with a virial mass M V - 1T ~ 2.9 x 10 15 /i _1 Mq (the most mas- 
sive object in the MS) and a merging (s ~ 0.1) cluster with 
M vir ~ 1.5 x 10 15 ft _1 M . The left panels show results for 
the GO simulation, the middle panels for the PC simulation 
and the right panels for the FO simulation. 

As has be en seen in previous simulations (e.g. 
iMotl et al.l 120051 ), the y distribution is very smooth. The 
most significant features are sharp edges associated with 
shocks; this is especially clear in the case of the merging 
cluster. Qualitatively, the maps look structurally similar be- 
tween models although their y values within a given pixel 
can be significantly different, with the GO and PC models 
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Figure 1. Top panels: Compton-?/ maps for the most massive cluster at z = in the MS (M v i r ~ 2.9 X 10 15 h —1 MQ, which we classify 
as regular) with X-ray surface brightness contours overlaid. Results are shown, from left to right, for the cluster in the GO, PC and 
FO models respectively. Images in the bottom panels are similar except a massive disturbed cluster (with M v ; r ~ 1.5 X 10 15 Ii~ 1 M.q) is 
shown. Each panel is 2r v ; r across and the scale and value of r v ; r is shown in the left-hand panels. The range of y values is given for each 
cluster in the scale at the bottom of each row; note the disturbed cluster has a lower maximum value than the regular cluster. The X-ray 
contours illustrate levels that are 10 per cent and 1 per cent of the maximum value in the map. The gross features are similar in all 3 
models for both clusters, although the X-ray maps reveal that the gas in the the PC clusters is the smoothest, while the GO clusters 
contain gas with the most small-scale structure. 



lying at either extremes. For the regular GO cluster, the 
mean y within the virial radius is (y) — 1.4 x 10~ 5 , with a 
range of values from 2 x 1CF 7 — 7 x 10~ 4 . For the PC clus- 
ter, the mean value is very similar although the maximum 
y (associated with the centre of the cluster) is almost half 
(ymax — 4 x 10 -4 ). This is due to the pre-heating of the gas 
that acts to smooth out the high density regions. 

It is also noticeable that the GO clusters contain a sig- 
nificant amount of small-scale structure in the gas. This is 
not clear in the y distribution but is evident from the over- 
laid X-ray surface brightness contours. These are clumps of 



low entropy gas associated with substructures in the cluster. 
Again, the pre-heating has smoothed these out by raising the 
entropy of the gas at high redshift. These features are also 
seen in the FO clusters, where heating is localised to haloes 
in which AGN feedback is occurring. 

2.6 Sky maps 

We also analyse simulated sky maps of the thermal SZ effect 
for the GO a nd PC models, using the stacked box approach 
pioneered by Ida Silva et al. (2000). This is an approximate 



3 X-ray surface brightness maps are calculated by replacing Xj 
in equation ( 1121 with piA(T;, Z), where pi is the density of hot 
gas particle i and Z = 0.3Z@ is the assumed metallicity. The 



cooling function, A(T, Z), is calculated for the soft [0.5-2] keV 
band. We normalise each surface brightness map to the maximum 
pixel value. 
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Figure 2. An example 5x5 square degree sky map of the Compton y parameter from the GO (left) and PC (right) simulations. The 
maps were smoothed with a Gaussian kernel with a full-width half-maximum of 1 arcmin. Circles illustrating virial radii of all clusters 
with M vir > 10 14 h -1 M s are overlaid. Although both maps contain the same clusters, the PC map clearly has a larger background and 
the y distribution within each cluster is smoother, as seen in the maps of individual objects. 



method for generating past light-cones using a finite num- 
ber of outputs. To do this we first compute the lookback 
time corresponding to a comoving distance of 50ft _1 Mpc. 
We then calculate successive lookback times, increasing the 
comoving distance in steps of A map = 100 h~ Mpc. These 
lookback times are used to find the nearest output time when 
simulation data are stored (a total of 160 snapshots were 
generated). We also calculate the comoving width required 
at each lookback time, corresponding to a fixed opening an- 
gle of # map = 5°. The final lookback time is chosen such 
that the comoving width is still smaller than the box-size, 
to avoid replication of the particles. The choice of 9 ml!Lp al- 
lows us to integrate the SZ effect out to a maximum redshift, 
Zmax = 4.7, using 47 snapshots; this is sufficiently large for 
the mean y signal to be converged in our simulations (see 
Fig- El discussed below). 

Once the required volumes are defined to make up the 
lightcone, the second stage is to use a random number gen- 
erator to construct a table of random translations, rotations 
(in steps of 7r/2 radians) and reflections about each of the 
three axes. This is done in order to minimise the chance of 
the lightcone containing the same cluster at different red- 
shifts (note the volume required at each time is always less 
than 20 per cent of the full simulation box because of our 
choice of A map ). The list of operations are then used to 
determine which particles are required to compute the con- 
tribution to the SZ signal from each redshift (used to create 
a so-called partial map) This stage is repeated 50 times to 
allow us to generate 50 quasi-independent realisations. 

The final stage is to generate the partial maps them- 



selves, by smoothing the appropriate gas particles on to a 
2D grid. This is done using the same technique as for indi- 
vidual clusters but now using a map area corresponding to 
#ma P x f? map at each redshift and a comoving depth of A map . 
Each partial map contains 1200 x 1200 pixels such that each 
pixel has an angular size, # p i x = 0.25 arcmin, comfortably 
smaller than the typical resolution of current SZ telescopes 
(1-10 arcmin). The 47 partial maps are then stacked for each 
realisation to make final maps of the y parameter. 

Fig. [2] shows an example Compton-y sky map for re- 
alisation 46, chosen because it contains a relatively large 
cluster. Both the GO (left panel) and PC (right panel) ver- 
sions are shown. The maps were smoothed using a Gaussian 
kernel with a full-width half-maximum of 1 arcmin, similar 
to the resolution of modern ground-based SZ telescopes such 
as SPT and ACT. 

The most striking difference between the two maps is 
the contrast: the PC map has a higher background than 
the GO map, making it harder to visually pick out the SZ 
sources associated with the clusters. This is due to the ex- 
tra thermal energy added to all the gas by the pre-heating 
process and can be quantified by measuring the mean y pa- 
rameter, averaged over all 50 realisations. For the GO run, 
we find (y) = 2.3 x 10~ 6 , increasing by more than a factor 
of four to (y) = 9.9 x 10~ 6 for the PC run. Although both 
values are below t he current constra int from COBi?/FIRAS, 
(y) < 1.5 x 10" 5 (|Fixsen et al.lll996T ). it is unlikely to be the 
case that the true background is as high as in the PC model, 
as this would erase many of the wea k neutral hydrogen ab- 
sorption lines seen towards quasars (|Theuns. Mo fc Schavd 
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Redshift 

Figure 3. Integrated contribution to the overall mean y signal 
from gas below a given redshift. The top panel is for the GO sim- 
ulation and the bottom panel for the PC simulation. Light curves 
are for individual maps and the dark curve is the average over 
all 50 maps. Each curve is normalised to the mean y averaged 
over all 50 maps, highlighting the scatter in the integrated signal 
between realisations. The simulations predict dramatically differ- 
ent redshift dependencies: the mean signal in the GO simulation 
comes from low redshfit (z < 2) whereas the opposite is the case 
for the PC simulation, due to the effect of pre-heating at 2 = 4. 



l200ll ; IShang, Crotts fc Haimar3l2007l : iBorgani fc Vielll200gh . 
The PC model therefore serves as an extreme test of the ef- 
fect of a high background although we will remove the mean 
y signal in our analysis in Section ^. 6l to mimic observations. 

The contribution to the mean y signal from gas at differ- 
ent redshifts is shown in Fig. The top panel shows results 
for all 50 maps in the GO simulation and the bottom panel 
for the PC simulation. Again, the difference between the 
two models is striking: the majority of the y signal comes 
from low redshift in the GO model (around 80 per cent 
from z < 2) whereas the opposite is true for the PC model 
(around 80 per cent from z < 3.5). Most of the mean y 
comes from overdense regions (groups and clusters) in the 
GO model that are more abundant at low redshift. In the 



PC case, most of the mean signal comes from m ildly over- 
dense gas at high redshift (|da Silva et al. 2001). Note also 
that the contribution from gas at z > 4 is approximately 
zero in the PC model, unlike in the GO case, where there is 
a small but non- negligible signal. This difference is due to 
the inclusion of radiative cooling in the former model which 
removes most of the (small amount of) ionised gas at these 
redshifts. 



3 HOT GAS PRESSURE PROFILES 

Fundamental to understanding the SZ effect from clusters 
is the hot gas pressure profile, since we can write the SZ Y 
parameter for a spherically symmetric cluster as 

F 60 o = — £j / P c (r)47rr 3 dlnr, (13) 

where P c = n c kT c is the electron pressure. The contribu- 
tion to Y500 will therefore be highest at the radius where 
r 3 P500 is maximal. If the gas is in hydrostatic equilibrium 
then the pressure profile ought to be structurally similar be- 
tween different clusters since it is directly constrained by 
the underlying gra vitational potential, which itsel f takes on 
a regular form (e.g. lNavarro, Frenk fc W hite 1997, hereafter 
NFW). 

We construct and compare spherically-averaged, hot 
gas mass- weighted pressure profiles using equation |TT}, for 
all clusters with M 500 > 1O 14 /i _1 M in our three (GO, 
PC and FO) models at z = 1 and z — 0. The profiles 
are re-scaled such that we plot dimensionless quantities 
x 3 P(x) / P500 against x, where x = r/rsoo and the scale 
pressure, P500 oc 

M Uo E ( z f /3 > is dete rmine d assuming a 
self-similar isothermal gas distribution l|Voitll2005l ). If clus- 
ters formed a self-similar population then these re-scaled 
profiles would be identical for both varying mass and red- 
shift. 

Median scaled profiles are shown in Fig. [4] split into low- 
mass (10 14 /i _1 M < M500 < 5x 10 14 /i~ 1 M ; triangles) and 
high-mass (M500 > 5 x 1O 14 /i _1 M0; squares) sub-samples. 
Comparing the high-mass clusters between the three models 
at z — 0, it is immediately apparent that the largest con- 
tribution to Y500 comes from radii close to rsoo, i.e. where 
P(r) oc r -3 . The profiles rise sharply (by around an order of 
magnitude) from the core outwards then stay level or gradu- 
ally decline at larger scales. The largest differences between 
the three models occur in the core region, where the PC and 
FO clusters have lower central pressures than the GO clus- 
ters due to the increase in core gas entropy from the extra 
heating. 

The low-mass clusters have very similar median profiles 
to the high-mass clusters in the GO s imulation, reflectin g 
the similarity of objects in that model l|Stanek et al.ll20lb1 ). 
In the PC and FO models however, the pressure profiles of 
the low-mass clusters have markedly different shapes from 
their high-mass counterparts. In particular, the scaled pres- 
sure in low-mass clusters is lower in the central region and 
is higher in the outer region, indicating that they are less 
concentrated than the high-mass clusters. Again this reflects 
the breaking of self-similarity caused by the feedback/pre- 
heating which has a larger effect in the lower mass clusters; 
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Figure 4. Scaled pressure profiles for clusters within the GO (top 
panels), PC (middle panels) and FO (bottom panels) simulations. 
Left panels are z = 1 results and right panels for 2 = 0. Me- 
dian profiles are shown separately for low-mass (10 14 h _1 M0 < 
M500 ^ 5 X lO 14 h _1 M0; triangles) and high-mass (M500 > 
5 X 1O 14 /i -1 M0; squares) respectively (note there are no high 
mass clusters at z = 1). The yellow (cyan) shaded band illus- 
trates the 16/84 per centiles (and thus represents the cluster-to- 
cluster scatter) for the high-mass (low-mass) sub-sample. Solid 
curves are best-fit generalised NFW profiles to the median pres- 
sure profiles. The vertical solid line represents the radius where 
gravity is softened in the cluster with the smallest rsoo (this is 
at a smaller radius than plotted for the FO model, for which a 
smaller softening was used, but we choose to use the same scale 
as in the other two models for ease of comparison). Note that the 
contribution to Y500 is predominantly from r > 0.5rsoo and so is 
not particularly susceptible to variations in the cluster core. 



the extra entropy given to the gas causes a re-distribution 
to take place, pushing the gas out to larger radius. 

Comparing the low-mass clusters at low and high red- 
shift, the GO model shows little evolution (the core pres- 
sures are slightly lower), while clusters in the PC model 
have significantly lower core pressures at z = 1. This re- 
flects the larger impact of the pre-heating on the gas at high 
redshift, since a cluster of fixed mass has a lower character- 
istic entropy at higher redshift from gravitational heating 
[K ~ M 2 / 3 E[z)~ 2 ' 3 ]. Interestingly, the scaled pressure pro- 
files in the FO model show little evolution with redshift, 
although the pressure in the outskirts (r > rsoo) is higher at 
2 = 0, reflecting the late-time heating of the gas by AGN. 

We also compare scaled pressure profiles between regu- 
lar (s 0.1) and disturbed (s > 0.1) clusters in Fig. O for 



Figure 5. As in Fig. [4]but the samples are now split into regular 
(s ^ 0.1; triangles) and disturbed (s > 0.1; squares) sub-samples. 
Disturbed (merging) clusters tend to have lower central pressures 
but higher peak values of x' i P(x). 



our samples with M500 > 10 14 /i _1 Mq. The largest differ- 
ences between the two sub-samples can be seen for the GO 
model, where the disturbed clusters (squares) have lower 
scaled pressure everywhere except around the maximum at 
r ~ 0.9r5oo- This is because the ongoing merger is com- 
pressing the gas (and therefore increasing its pressure) at 
large radius while the inner region has yet to respond to the 
increase in the mass of the system. Note that since Y500 is 
proportional to the area under the pressure profile, there will 
be a noticeable offset in the Y500 — Mmo relation, where a dis- 
turbed cluster has a smaller Y500 than a regular cluster with 
the same mass (see the next section). These differences are 
still present but at a lower level in the PC and FO models, 
where the higher entropy of the gas in lower-mass clusters 
means that it is less easily compressed. This in turn leads to 
a negligible offset between regular and disturbed clusters in 
the Y500 — M500 relation, as we will show in the next section. 

The shaded bands in Figs. [4] and \5\ illustrate the 16/84 
per centiles for the two respctive sub-samples and thus gives 
an indication of the cluster-to-cluster scatter. We show this 
more clearly in Fig. [6] where we have normalised the clusters 
in the low and high mass sub-samples to the generalised 
NFW model that best fits the median profile (see below). 
Although the scatter at fixed radius is quite low compared 
with some other properties such as X-ray surface brightness, 
it is nevertheless appreciable and can be as high as 30-50 per 
cent beyond rsoo • Thus it is clearly not accurate to assume a 
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Figure 6. As in Fig. [4] but the median pressure profiles (and 
scatter) are now shown relative to their best-fit generalised NFW 
model, clearly showing the size of cluster-to-cluster variations 
(that can be as large as 50 per cent). Solid and dashed curves 
are observed mean pressure pro files from low-redshift X-ray data 
(REXCESS: lArnaud et alj20 1Q). again scaled to our best-fit gen- 
eralised NFW model profiles, assuming the median mass from our 
low and high-mass sub-samples respectively. In the outer regions 
(r > 0.5r5oo) the high-mass clusters in the PC and FO models 
fit the latter profile quite well (to within 10 per cent or so), but 
the difference is larger for low-mass clusters, especially in the core 
regions. 



single profile to describe all clusters, especially around rsoo 
and beyond, where much of the SZ signal comes from. 



3.1 Generalised NFW model 

In a previous st udy of hot gas pressure profiles in cos molog- 
ical simulations, iNagai. Kravtsov fc Vikhlininl (|2007t ) found 
that the mean pressure profile of their simulated clusters 
could be well described by a generalised NFW (GNFW) 
model with five free parameters 



P(r) 







(14) 



p500 ^(l+U^W-T)/ ' 

where u = csoo^, C500 is the concentration parameter, Po 
the normalisation parameter and (7, a, /3) determine the 
shape of the profile at small [u <C 1), intermediate (u ~ 1) 
and large (u 3> 1) radius respectively. The GNFW model 
has been shown to provide a good descriptio n to the pres- 
sure p rofiles of X-ray g roups and clusters (e.g. lArnaud et all 
l201fj ; ISun et all l201lT ) and is being used to optimise SZ 



Table 2. Best-fit parameters for the generalised NFW model 
when applied to our median hot gas pressure profiles. Column 1 
gives the redshift; column 2 the simulation model and cluster sub- 
sample (LM and HM refer to the low and high-mass sub-samples, 
respectively); and columns 3-7 the parameter values (see text for 
their meanings). 



Redshift 


Clusters 


Po 


£500 


7 


a 


/8 


2 = 


GO/LM 
PC/LM 
FO/LM 


33.788 
6.317 
4.732 


2.925 
0.517 
1.052 


0.267 
0.090 
0.298 


0.944 
0.901 
1.108 


1.970 
1.603 
2.371 




GO/HM 
PC/HM 
FO/HM 


6.756 
0.938 
3.210 


1.816 
0.183 
1.974 


0.519 
0.584 
0.605 


1.300 
1.114 
2.041 


2.870 
11.885 
2.989 


2 = 1 


GO/LM 
PC/LM 
FO/LM 


11.994 
0.856 
2.734 


0.700 
0.539 
0.349 


0.345 
0.512 
0.375 


0.837 
1.447 
1.055 


3.610 
4.038 
5.049 



cluster detection in data f rom the Planck satellite (e.g. 
iPlanck Collaboration! |2011ch . 

We have applied the GNFW model to our simulated 
clusters and the results for the median profiles can be seen as 
solid curves in Figs. [4] and [5] We also normalise our pressure 
profiles to the best-fitting median GNFW profile in Fig. [6] 
The residual values for our median profiles are also shown 
(as triangles and squares for our low and high mass sub- 
samples) and are clearly at the per cent level. Such small 
residuals are not surprising given the model contains five free 
parameters (once rsoo is specified). The best-fit parameter 
values themselves are listed in Table [2] 

To investigate the distribution of GNFW parameters 
and any degeneracies that arise between parameters, we plot 
marginalised likelihood distributions for the FO model at 
2 = in Fig. [7] The full five-dimensional likelihood dis- 
tribution is estimated by fitting the GNFW model to in- 
dividual clusters and computing the frequency of parame- 
ters [log 10 Po, log 10 C500, 7, ct, 0\ over a five-dimensional grid, 
which is then normalised such that the sum over all allowed 
parameter values is unity. We assume, as prior information, 
that the allowed range for each parameter is as specified on 
the axes in Fig. [7] and that each value is equally likely. 

The diagonal panels in Fig. [7] show the marginalised 
ID likelihood distributions for each of the five parameters, 
while the off-diagonal panels show the 68 per cent confi- 
dence regions for the full range of marginalised 2D distribu- 
tions, smoothed to reduce noise. The concentration param- 
eter is strongly correlated with the normalisation parameter 
but does not correlate strongly to any of the slope parame- 
ters. Interestingly, the normalisation is anti-correlated (and 
therefore degenerate) with the slope parameters. Finally, the 
slope parameters are strongly correlated with one another. 
It is therefore clear that a simpler model with fewer slope 
parameters could potentially be found that describe these 
simulated data. However, the flexibility of the GNFW model 
allows a wide range of profiles to be accurately described us- 
ing a simple formula. This is especially true when cool-core 
clusters are included; these are absent in our current models 
and so we plan to extend our work to investigate cooling 
effects in a future study. 
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Figure 7. Estimated likelihood distributions for the best- fit generalised NFW model parameters that describe the hot gas pressure profiles 
in the FO simulation at z = 0. The off-diagonal panels show 68 per cent confidence regions for the marginalised 2D distributions, for all 
parameter combinations. The diagonal panels show the ID marginalised distributions for each parameter, with the best-fit parameter 
values from REXCESS data overlaid as dashed vertical lines. The three slope parameters (a,/3,7) are strongly correlated with one 
another and are all negatively correlated with the normalisation, log 10 Po- The concentration parameter, log 10 C5oo, is only correlated 
with the normalisation. 



3.2 Comparison with observational data 

We also compare th e simulated pr o files w ith the pressure 
profile presented by lArnaud et alj l|2010t ). compiled from 
low-redshift X-ray observations (for r < rsoo; the REXCESS 
sample) and other numerical simulations (for r > rsoo). It 
therefore provides information on the realism of our simu- 
lated pressure profiles as well providing a useful comparison 
with other simulations (on large scales). 

The lArnaud et all (|201Ch profile is based on the GNFW 
model, modified to account for additional (weak) mass de- 
pendence in the observational data 



P(x) = Pgnfw 



M 5 , 



3 x 10 14 M Q 



(15) 



where Pgnfw is the GNFW pressure profile given 
in equation (|14[l with parameters, [Po, C500, 7, ex, P] = 
[8.403, 1.177, 0.3081, 1.0510, 5.4905] and a P = 0.12. We show 
this profile, evaluated for the median mass values of our two 
sub-samples, in Fig. [S] the dashed curve is for the low-mass 
sample, plotted relative to our best-fit GNFW profile, while 
the solid curve is for the high-mass sample. The Arnaud et 
al. parameters are also shown as dashed lines in Fig. E] 

Comparing with our z — results, as is most appropri- 
ate, the median GO profiles agree to within 30 per cent or so, 
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over the plotted range of radii and for both mass ranges. For 
the PC and FO clusters, the agreement is very good at large 
radius (> 0.5rsoo) for high-mass clusters, where the Arnaud 
et al. profile is only around fO per cent higher and within 
the intrinsic scatter of our simulated profiles. The low-mass 
clusters are more discrepant, with the steeper Arnaud et al. 
profile being 20-30 per cent lower at rsoo- This suggests that 
our simulated clusters contain gas that is at higher pressure 
at rsoo than in those used for the Arnaud et al. profile at 
large radius. Given that the feedback in our models is likely 
to be stronger than in the simulations used in the Arnaud 
et al. study, this discrepancy in pressure is probably due to 
the effects of radiative cooling, absent in our models and 
likely significant in the other simulations (see the discus- 
sion in Section f4. 41 we also note that Arnaud et al. already 
corrected for the effects of baryon fraction). Even larger dif- 
ferences are present in the inner regions; there, the Arnaud 
et al. profile is significantly higher than our simulated re- 
sults. Again, cooling is the likely culprit here as its effect is 
strongest in the densest regions. 

An important uncertainty in the observed profile esti- 
mation is the effect of hydrostatic bias, i.e. systematic offsets 
in rsoo and M500 from their true values, when estimated from 
the equation of hydrostatic equilibrium. As we will show in 
Section [5l hydrostatic mass is biased low with respect to 
the true mass and is most significant for the GO model (the 
estimated-to true mass ratio is around 0.7 for the GO model, 
compared with around 0.9 for the PC/FO models). The ef- 
fect of this bias is to increase the scaled pressure at fixed 
scaled radius, as both the scale radius, rsoo, and the scale 
pressure, P500 oc decrease, on average. We discuss the 

effect of hydrostatic bias on the Y500 — A/500 relation in detail 
in Section but note here that we have explicitly checked 
how this affects the pressure profiles for each model. To do 
this, we first re-defined our sub-samples using the estimated 
masses. We then compared the shift in pressure at the es- 
timated value of rsoo between the median scale profile and 
the Arnaud et al. profile, for both low-mass and high-mass 
sub-samples. We also re-computed the pressure profiles us- 
ing the spectroscopic-like temperature, rather than the hot 
gas mass-weighted temperature, as this will be closer to the 
X-ray temperature profile used by Arnaud et al. 

We find that the combined effect of these changes is 
largest for the GO model, where the median pressure pro- 
files from both sub-samples are now within fO per cent of 
the Arnaud et al. values at rsoo- The increase in the scaled 
pressure profile due to hydrostatic bias is counteracted by 
a decrease due to the use of spectroscopic-like temperature, 
which is lower than the mass-weighted temperature for this 
model (see Section f4.2p . The two effects are smaller for the 
PC and FO models and so we see very similar results to 
those before these changes were applied. Thus, the scaled 
pressure profiles for the low-mass clusters in these models 
are still around 30 per cent lower than the Arnaud et al. 
profile at rsoo. 



suits with recent simulations before going on to consider the 
effect of projection of large-scale structure along the line-of- 
sight. The effect of hydrostatic bias on the scaling relations 
will be considered in the next section. 

We consider the scaling relations between Y500 and sev- 
eral other properties: the total mass, M500; the hot gas mass, 
A/gas, 500; the X-ray spectroscopic-like temperature, T B \; and 
the analagous X-ray quantity to Y500, Yx,5oo- We note that 
the Yx — M500 rel ation (not conside red here) has already 
been presented bv lShort et al.l (|201 Ol ) and scaling relations 
for the lo wer density contrast, A = 200, fo r the GO and PC 
models by IStanekT Rudd fc Evrardl (|2009l ). Fl 

We follow the standard procedure and assume that the 
mean relationship between Y500 and the independent vari- 
able can be adequately described by a power law and is thus 
a linear relationship in log-space. □ We estimate the slope 
and normalisation of the relation by performing a least- 
squares fit to the data 



E(zyY 500 = 10 A (X/X ) B . 



(16) 



where A and B describe the best-fit normalisation and slope 
respectively and Xo is the pivot point, suitably chosen to 
minimise co-variance between the two parameters. For the 
power-law index 7 we choose the appropriate value for self- 
similar evolution, so if our clusters evolve self-similar ly we 
should see no change in the best-fit parameters A and B. 
We also estimate the scatter in log 10 (Y 5 oo), vio Sl0 y, as 



0-io glo y 



\ 



1 



N ■ 



^2 [logio Yi - log 10 Y hf (Xi 



(17) 



where the index i runs over all N clusters included in the 
fit and Ybf is the best-fit Y500 value for a cluster with prop- 
erty, Xi. Note that the scatter in \nY is simply eri n y = 
ln(f0) crio gl0 y. 



4.1 The Y500 — M500 relation 

The most important scaling relation is that between SZ 
flux and mass. We present our Y500 — M500 relations in 
Fig. [5] for the GO model (top panels), PC model (mid- 
dle panels) and FO model (bottom panels). Results are 
shown both for z — 1 (left panels) and z — (right pan- 
els). The best-fit relation to all clusters in each panel with 
fO 14 /i _1 M < M 50 o < 10 15 /i _1 M Q is shown as a solid line 
(best-fit parameter values and their uncertainties are listed 
in Table 13). 

ft is clear that there is a very tight correlation between 
Y500 and M500 in all three models at both low and high 
redshift. At z = the intrinsic scatter about the best-fit 
power law relation is only ~ 4 per cent, with sub-per cent 
variations between models, making this particular relation 
one of the tightest known cluster scaling relations involving 



4 SZ SCALING RELATIONS 

We now present SZ scaling relations for our simulations and 
compare them specifically with the recent analysis of data 
from Planck and XMM-Newton. We will also compare our re- 



4 We have independently verified that our GO and PC results, 
when using A = 200, are c o nsiste nt with theirs, but as was 
pointed out by IViana et al. (2011) the Y200 values given in 
IStanek. Rudd fc Evrardl d2009h quoted with incorrect units. 

5 IStanek et al. present quadratic fits to the PC data but 
we find this only to be important when the lower-mass groups are 
included, as was the case in that study. 
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Table 3. Best-fit parameters for simulated SZ scaling relations at z = and 2 = 1. Column 1 gives the scaling relation being considered; 
column 2 the pivot point (in appropriate units); column 3 the redshift; column 4 the simulation model; and columns 5-7 the best-fit 
values for the normalisation, power-law index and scatter in log 10 Y500 respectively. Quoted uncertainties correspond to either the 16th 
or 84th per centile (whichever is largest), estimated using the bootstrap re-sampling technique. 



Relation 








Kcdsmtt 


Model 


* 

A 


B 




E(z)- 2 / 3 Y 500 


- M 500 


3 x 10 14 h~ 




z = 


GO 
PC 
FO 


-4.754 ±0.002 
-4.774 ± 0.003 
-4.744 ± 0.003 


1.670 ± 0.007 
1.794 ± 0.009 
1.69 ± 0.02 


0.041 ±0.001 
0.045 ± 0.001 
0.043 ± 0.003 










z = l 


GO 
PC 
FO 


-4.79 ±0.01 
-4.82 ±0.01 
—4.75 ± 0.01 


1.60 ±0.04 
1.83 ± 0.04 
1.63 ± 0.04 


0.048 ± 0.003 
0.059 ± 0.005 
0.037 ± 0.004 


E{z)- 2 t 3 Y M0 


— M gaS! 500 


3 x 10 13 h~ 




2 = 


GO 
PC 
FO 


-5.098 ±0.001 
-4.887 ± 0.001 
-4.889 ± 0.003 


1.650 ±0.007 
1.478 ± 0.008 
1.45 ± 0.01 


0.029 ±0.001 
0.018 ±0.001 
0.025 ±0.002 










Z = 1 


GO 
PC 
FO 


-5.145 ±0.004 
-4.844 ± 0.006 

rr nnT i n nn a 

— 5.007 ± 0.004 


1.61 ±0.03 
1.46 ± 0.05 

1 q 1 n no 

1.53 ± 0.03 


0.034 ±0.003 
0.016 ±0.005 

n no 1 n one 

0.028 ± 0.006 


E(z)Y 500 - T s 


1 


5keV 




z = 


GO 
PC 
FO 


-4.27 ±0.03 
-4.706 ± 0.006 
-4.665 ± 0.006 


2.5 ± 0.2 
3.16 ± 0.04 
3.11 ± 0.07 


0.19 ±0.02 
0.060 ± 0.002 
0.078 ± 0.007 










z = l 


GO 
PC 
FO 


-4.03 ±0.07 
-4.910 ±0.006 
-4.54 ±0.02 


3.0 ±0.6 
3.38 ± 0.06 
2.8 ± 0.1 


0.16 ±0.07 
0.047 ±0.004 
0.086 ± 0.009 


E(z)- 2 / 3 Y 50 o 


- E(z)- 2 / 3 Y X:500 


1 x 1CT 5 h' 


" 2 Mpc 2 


z = 


GO 
PC 
FO 


-4.952 ±0.003 
-5.020 ± 0.001 
-5.012 ±0.001 


1.049 ± 0.008 
1.002 ± 0.002 
0.998 ± 0.004 


0.058 ± 0.002 
0.015 ±0.001 
0.018 ±0.001 










z = 1 


GO 
PC 
FO 


-4.882 ± 0.009 
-5.015 ±0.002 
-5.007 ± 0.003 


1.05 ± 0.03 
0.999 ± 0.004 
0.99 ± 0.01 


0.052 ± 0.004 
0.009 ± 0.001 
0.020 ± 0.004 



gas; this finding is c onsistent with previous simulation s with 
fewer clusters (e.g. Ida Silva et alJ liooj ; iNagail l200rj ). The 
distribution of residual Y500 values about the best-fit relation 
is well described by a log-normal distribution of width a — 
crinY (Fig.[9t. This is in agreement wi th previous work (e.g. 
IStanek et ai1l2010l ; iFabian et a"ill201ll ). 

The normalisation of the 2 = relation also varies very 
little between models, the maximum variation being around 
7 per cent. The best-fit slope also varies by around 7 per 
cent, from 1.67 in the FO model (very close to the self- 
similar value of 5/3) to 1.79 in the PC model. As discussed in 
IShort fc Thomas! (|2009l ) for the Fx - M 500 relation, the sim- 
ilarity between the models can be explained by the increase 
in gas temperature compensating for the drop in gas mass, 
required to maintain virial equilibrium (since Y cx Af gas T 
and is thus proportional to the total thermal energy of the 
gas). The agreement between the GO and PC/FO models 
is better here than for the Yx — A/500 relation as Yx is de- 
fined using the spectroscopic- like temperature, T s \, which is 
weighted more heavily by low entropy gas; we discuss this 
point further below. 

We have also investigated the dependence of the Y500 — 
M500 relation on redshift. In the left-hand panels of Fig. [HI 
we present results for z = 1, allowing a simple comparison 
to be made with the 2 = results for each model. It is 
evident that the clusters evolve close to the self-similar ex- 
pectation in all three models, given that the normalisation 
and slope changes very little between the two redshifts (see 
also Table [3|. To quantify this further, we have also plotted 



the best-fit normalisation, slope and scatter as a function 
of redshift in Fig. 1101 where we used all available outputs 
from 2 = to 2 = 1. (Equivalent plots for the other scaling 
relations are provided in the Appendix.) 

The dependence of the best-fit slope with redshift for all 
three models is shown in the top panels of Fig. llUI For clarity 
we normalise the slope to the median value at 2 < 0.3 and 
the yellow bands indicate the uncertainties (using the 16/84 
per centile values). All three models are consistent with no 
evolution in slope to z = 0.3, then some mild evolution is 
seen at higher redshift, where the number of massive clusters 
drops. This evolution is very minor, however, as the slope 
remains within around 5 per cent of the low redshift value. 

The variation in normalisation with redshift is shown in 
the middle panels of Fig. [TU] Here, we have fixed the slope 
at the 2 < 0.3 median value and just allowed the single nor- 
malisation parameter to vary. Again, we factored out the 
self-similar evolution and normalised to the 2 = result, 
so a value consistent with zero corresponds to self-similar 
evolution. In the GO and PC cases, the normalisation is 
consistent with self-similar evolution to 2 = 0.3, afterwards 
there is some negative evolution (i.e. the relation evolves 
slightly more slowly than predicted from the self-similar 
model), especially in the PC case. The FO model shows 
different behaviour: at low redshift (2 < 0.3), Y500 increases 
more rapidly with redshift than the self-similar case (at fixed 
mass), then at higher redshifts evolves in accordance with 
the self-similar expectation. These differences in evolution 
are likely to be real and reflect the varying gas physics. In 
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Figure 8. Scaling relations between the SZ flux, Y500, and to- 
tal mass, M500, for the GO (top panels), PC (middle) and FO 
(bottom) models at z = 1 (left) and z = (right panels). The 
Y500 values at z = 1 are re-scaled such that no change in the 
relation corresponds to self-similar evolution (Y500 <x E(z) 2 / 3 
at fixed mass). The solid diagonal line is a least-squares fit to 
the relation. The best-fit power-law to z < 0.5 Planck/ XMM- 
Newton data jPlanck Collaborationboildl is shown in all panels 
as a dashed line, while the box illustrates the intrinsic scatter in 
the observed relation. Triangles represent disturbed clusters (with 
s > 0.1) while squares are regular clusters. The results are very 
similar for all three models and there is no evidence for significant 
departure from self-similar evolution. 



the GO case, the gas at high redshift is slightly colder than 
expected (due to an increase in the merger/accretion rate 
leading to a larger residual unthermalised component). In 
the PC case, the high redshift pre-heating leads to a deficit 
in gas mass but the clusters start to recover at low red- 
shift as the entropy scale at fixed mass set by gravitation 
is larger. Finally, in the FO case, the feedback from black 
holes is stronger at late ti mes, leading to a decrease in gas 
content (|Short et al.ll201ol ). In all three cases, however, the 
effect on the normalisation is still small; the largest change 
is from the PC model at z = 1 where only a 10 per cent 
decrease is seen. 

Finally, we illustrate how the scatter in the Y500 — A^soo 
relation evolves with redshift in the bottom panels of Fig. 1101 
The z = value is also shown as a dashed horizontal line for 
clarity. Again, the picture is consistent with minimal change; 
the scatter only increases to z — 1 by 0.01 or so in the GO 
and PC cases, and decreases by less than 0.01 in the FO 
case. 



Figure 9. Distribution of residual Y500 values about the best-fit 
Y500 — A^500 relation, plotted using natural logarithms, for each 
of the models at z = 1 and z = 0. A normal distribution of 
width cr = a\ n y is overlaid; it is clear that this provides a good 
description of the scatter. 



4.2 Relationship between Y500 and observables 

We also present scaling relations between Y500 and other 
key X-ray observables. Fig. [11] shows Y500 — M gas ,500 rela- 
tions, laid out as before. This relation is interesting to study 
because it essentially probes non self-similar behaviour in 
the mass- weighted temperature, T m , of the gas, since Y oc 
AfgasTm and thus M gas appears on both axes. Here we 
fit data within the range, 2 x 10 13 h~ 1 M Q < M gaa , 500 < 
2 x 10 14 h, _1 M . As with the Y500 — M500 relation, the slope 
from the GO model at z — is close to the self-similar value 
of 5/3. The PC and FO models have shallower slopes, due 
to the increase in the temperature of the gas in low-mass 
clusters. As might be expected, the scatter in the relation is 
even tighter than for the Y500 — Aftoo relation, and is typi- 
cally 0.02-0.03. The distribution of the scatter is also close 
to log-normal. From comparing the z — 1 and z — results, 
both GO and PC models predict evolution that is close to 
self-similar (the normalisation is within 5 per cent of the 
z — value out to 2 = 1) but the FO relation evolves more 
slowly with redshift (~ 10 per cent lower at 2 = 1), again 
due to the increase in feedback from the AGN at late times 
that additionally heats the gas. This evolutionary behaviour 
is confirmed when studying the relation at all available red- 
shifts from z = to z = 1, in Fig. IA1I which also shows that 
the slope and scatter vary little. 

We also consider scaling relations between Y500 and X- 
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Figure 10. The dependence of the slope, normalisation and scatter of the E(z)~ 2 / 3 Y$oo — M500 relation with redshift, for the GO 
(left), PC (middle) and FO (right) simulations. Results are plotted from z = to z = 1. In the top panels, the best-fit slope at each 
redshift is normalised to the median slope for outputs at z < 0.3 (shown as a vertical line). The middle panels illustrate the redshift 
dependence of the normalisation after the self-similar dependence has been taken out; the normalisation is divided by the 2 = value in 
this case. In the bottom panels, the rms scatter in log 10 Y, o"i ogl0 y-, is shown as a function of redshift. For both the normalisation and 
scatter values, the slope was fixed to the z < 0.3 median value when performing the fits. The bands in all panels illustrate 16 and 84 
per centiles, calculated by bootstrap resampling the data. All three models predict very little evolution in the slope and normalisation 
of the E(z)~ 2 / 3 l5oo — -^500 relation to z = 1 and the intrinsic scatter remains small (c logl0 y < 0.06). 



ray spectroscopic-like temperature, T s i, and show results in 
Fig. 1121 with the redshift dependence of the slope, normali- 
sation and scatter illustrated in Fig. IA2I Here, we further re- 
strict our sample to contain only clusters with kT s \ > 3 keV, 
as the spectroscopic-like temperature only applies to hot 
clusters where thermal bremsstrahlung dominates the X-ray 
emission. This reduces our samples to 136 (12), 583 (102) 
and 179 (73) clusters at z = (z = 1) in the GO, PC and 
FO models respectively. Note the more severe reduction in 
the GO case; the non-gravitational heating in the PC and 
FO models increases T B \ at fixed mass, relative to the GO 
case, and thus increases the number of clusters in their X- 
ray temperature-limited samples. Best-fit relations are then 
calculated for clusters in the range, 3keV < kT B \ < lOkeV. 

The GO model relation has a slope that is consistent 
with the self-similar expectation (B = 5/2) at z = and 
2 = 1. The relation evolves slightly faster than the self- 



similar model (the normalisation is around 10 per cent 
higher than expected at z = 1), while the scatter is approxi- 
mately constant at all redshifts, but is much higher than for 
the previous relations (o"i oglo y ~ 0.15 — 0.2). This last point 
is due to Tsi being a much noiser property as it is sensitive to 
the clumpy, low entropy gas that is prevalent in this model. 
We also note that the scatter is poorly described by a log- 
normal distribution. In comparison, the PC and FO models, 
which have much smoother gas, typically have lower scatter, 
o"io glo y — 0.05 — 0.1, that is well described by a log-normal 
distribution. The slope in these two models is significantly 
steeper (B ~ 3) and the evolution of this relation shows 
the largest departure from self-similarity (up to 20 per cent 
lower/higher at z = 1 in the PC/FO models). 

Finally, in Fig. [T3] we plot Y500 against Fx, 500 for our 
cluster samples and show the redshift dependence of the 
scaling relation parameters in Fig. IA3I We do this to di- 
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Figure 11. Scaling relations between Y500 and hot gas mass, 
Afgas,500 for the GO, PC and FO models at 2 = 1 and z = 
0. Details for each panel are the same as described in Fig. \E\ 
The scatter is particularly small for the simulated version of this 
relation as the quantities plotted on the two axes are strongly 
correlated. 



Figure 12. Scaling relations between Y500 and X-ray 
spectroscopic-like temperature, T s i, evaluated outside the core. 
Details for each panel are the same as described in Fig. [8] This 
relation has relatively large scatter for both the observations and 
simulations, since T s ] is a noisier property than the other observ- 
ables being considered in this study. 



rectly highlight how the choice of gas temperature affects 
the results: any deviation from Vsoo = Yx,500 must be due to 
differences between mass and X-ray weighted temperatures. 
No significant deviation is seen in the PC and FO mod- 
els (the difference in normalisations at E(z)~ 2 ^ 3 Yx,soo = 
1CF 5 YT 2 Mpc 2 is less than 5 per cent) and there is very lit- 
tle scatter (oi og y — 0.01 — 0.02) at low and high redshift, 
that again has a distribution that is log-normal. The GO 
model, on the other hand, shows a significant bias, such that 
Y500 — lTYx.soo at z = 0, increasing to Y500 — l-31x,500 at 
z — 1. The scatter is also significantly larger than for the 
other two models, o"i ogl0 y — 0.05, and the distribution is 
skewed to lower values. Again, these results demonstrate 
that the clumpier gas in the GO model has a stronger effect 
on the X-ray properties than the SZ properties. As we shall 
see in Section [5] this has important consequences for our 
hydrostatic mass estimates. 

4.3 Effect of dynamical state 

It is also interesting to consider whether clusters undergo- 
ing mergers are offset from the main Y500 scaling relations 
as they could add to the intrinsic scatter. We mark our dis- 
turbed (s > 0.1) sub-samples as triangles in each of the fig- 
ures presenting scaling relations, discussed above (Figs. [8} 



I13p . Note that while a large value of s is indicative of an 
ongoing merger, not all dynamically disturbed c lusters have 
large values of s (jRowlev. Thomas fc Kavll2004 ). 

As predicted from studying the hot gas pressure profiles 
in Section [3] the only significant offset seen between regular 
and disturbed clusters is for the GO model, where disturbed 
objects lie slightly below the Y500 - M 500 and Y500 — M gaSi5 oo 
relations, and above the Y500 — Yx.soo relation (there are 
not enough disturbed clusters to say anything conclusive 
for the Y500 — T s \ relation). This suggests that there is a 
significant difference in the fraction of unthermalised energy 
between regular and disturbed clusters in this model. In the 
case of the Y500 — M500 and Y500 — M gas ,5oo relations, the 
mass- weighted temperature is lower for disturbed clusters of 
the same mass than regular clusters, leading to the negative 
offset. The effect is exacerbated when T s \ is considered (since 
it is weighted towards the cooler component), leading to a 
positive offset in the Y 5 oo — Yx,soo relation. 

4.4 Comparison of Y500 — M500 relation from other 
simulations 

Given the importance of the Y500 — M500 relation for cosmo- 
logical applications and its apparent insensitivity to cluster 
gas physics, it is important to compare our results to those 
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Figure 13. Scaling relations between the SZ Y500 and the X-ray 
analogue, lx,500- Details for each panel are the same as described 
in Fig. [8] The dot-dashed line corresponds to Y500 = ^x,soo; only 
the GO simulation shows a significant offset from this relation, 
due to the presence of clumpy low-entropy gas. 



from other groups using different simulations. A number of 
studies have been performed with varying assumptions for 
both the cosmology and gas physics, as well as the use of dif- 
ferent algorithms for the Af-body and hydro d ynamics solvers 
(e.g. IWhite. Hernquist fe Springell I200I Ida Silva et all 



20041; llvlotl et alj 120051; INagail l200rj; iBonaldi et all 



Hallman et al 



2007; 



2009 



Aghanim, da Silva fe Nunes 
Yang. Bhattacharva fe Rickerl 120101 ; iKrause et al, 2011 
Battaglia et al.ll201ll ). 



2007 



W e choo se to compare our r esults with the work of 
INagail (|2006t ) and IKrause et al] (|201lh for two reasons. 
Firstly, both groups presented results for A = 500 and are 
thus most readily comparable with ours. Secondly, the two 
groups used very different codes, so it is useful to also include 
that uncertainty in our comparison. In Fig. [TJ we compare 
our best-fit Y500 — M500 relation at z — from the FO model 
(solid line) with the results of these authors. To highlight the 
differences between simulations, we normalise all results to 
our best-fit FO relation. We also have to make a correction 
for the different baryon fractions used in the simulations, 
since Y oc /b. In both cases, the bar yon fraction is lower 
than our ado pted value of fb = 0T8 (|Nagaill200rj assumed 
/ b = 0.14 and lKrause et alfclll assumed f b = 0.133). Note 



that this is not a perfect correction as it does not account 
for the non-self-similar behaviour of the baryon fraction with 
cluster mass. 
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Figure 14. Comparison between the FO Y500 — A^500 scaling 
relation and results from other simulations at z = 0. All results 
are normalised to the best-fit FO I500 — M500 relation and Y500 
values from other simulations have been re-scaled to account for 
differences in baryon fraction (see text for details). The square 
symbols are individual cluster Y500 values from the FO model and 
the dashed line the best-fit GO relation. The dot-dashed lines are 
best-fit relations from INagail <2006ft : the upper line corresponds 
to their non-radiative (AD) simulation and the lower line their 
run with cooling and star for mation (CSF). Finall y, the triple- 
dot-dashed lines are taken from IKrause et ah I <201ll) for their re- 
stricted A and B samples (upper and lower lines respectively). 



INagail (|2006t ) presented results for 11 clusters simu- 
lated with the art code (e.g. iKravtsov. Klypin fc Hoffman! 
2002), that uses the adaptive mesh refinement technique to 
model hydrodynamics. Two sets of runs were studied, a non- 
radiative run (labelled AD) and a run with cooling and star 
formation (labelled CSF). Out of the 11 clusters, 6 have 
M500 > 10 14 /i _1 M (c.f. our FO model with 188 clusters 
in this mass range). The upper dot-dashed line in Fig. 1141 is 
their best-fit relation to the AD clusters. The slope of their 
relation (1.66T0.09) is in agreement with our (non-radiative) 
GO result (1.670 ± 0.007; dashed line) while the normalisa- 
tion is within 10 per cent of ours. Such good agreement is 
reassuring given the different hydrodynamic methods used, 
although the large difference in sample size must be borne 
in mind. Their CSF result is shown as the lower dot-dashed 
curve in Fig. 1141 comparing with our FO relation, their nor- 
malisation is significantly lower (20-30 per cent). As the au- 
thor points out, the reduction in SZ signal in the CSF run is 
mainly due to the lower gas fraction caused by (over-)cooling 
and star formation that removes hot gas from the ICM. As 
we discussed earlier, the gas fractions in the FO run are also 
lower than in the non-radiative case but the mechanism re- 
sponsible (strong feedback) compensates for this by heating 
the gas to a higher temp erature. 

IKrause et alj l|201ll ) present results for two cluster sam- 
ples, A and B, shown as the upper and lower triple-dot- 
dashed lines in the figure. Both samples were simulated with 
the same GADGET2 A r -body/SPH code as used in this study 
but contained different assumptions for the gas physics. 
Sample A contained 39 clusters re-simulated from a large 
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parent volume while sample B was a mass-limited sample of 
117 objects, taken fro m a sing l e simu lation. While both sam- 
ples are larger than in iNagail (|2006h the number of massive 
clusters is still significantly smaller than in our FO sample. 
The two samples (we show results restricted to clusters with 
M500 > 2 x 10 14 ft _1 M0) compare well with ours once the 
different baryon fraction is scaled out. The normalisation in 
both cases is within 10 per cent or so, although the slope is 
slightly flatter, a result that appears only marginally signif- 
icant (the slope of sample B is 1.64 ± 0.03 compared with 
the FO slope of 1.69 ± 0.02). 



4.5 Comparison with observational data 

We have also compared our results to observational data, 
now that blind SZ surveys are starting to yield sig- 
nificant numbers of (SZ-selected) clusters, enabling es- 
timates of the Ycm — -M500 relation to be per formed 
|Andersson et al.ll201l] ; iPlanck Collaborationll201ld 'l. Here, 
we compar e our results with those fro m the Planck Col- 
laboration l|Planck Collaboration! 1201 id . hereafter PXMM), 
although we note that their best-fit Y500 — M500 relation is 
similar to the SPT result derived from a lower number of 
clusters bv lAndersson et al.l (|20 1 lh . 

The PXMM sample consists of 62 clusters with z < 
0.5 and used X-ray data from XMM-Newton to define the 
size (rsoo) and mass (M500) of each cluster, calibrated us- 
ing Jhe_X^ray_M5qo_j-Yx, 500 relation previously derived 
by lArnaud et all (|2010h . Once the cluster size was de- 
fined, the SZ flux was measured using a multi-frequency 
matc hed-filter technique, based on the ICM pressure pro- 
file o f lArnaud et all (|2010h . We show their best-fit results to 
the Y500 - M 50 o, Y 500 - M gas ,5oo and F 50 o - T s i relations as 
dashed lines and illustrate their intrinsic scatter with boxes, 
in Figs. 1 81121 (Note we show these in both panels to help 
gauge the sense of evolution in our simulated relations, but 
the observed fits are more applicable to our z = results.) 

It is remarkable how well the PXMM results agree with 
our PC and FO models; only the Y500— T s \ relation shows any 
obvious offset but that is nevertheless small. The reason for 
such good agreement is not obvious or necessarily expected, 
given the complicated procedure involved in deriving the 
observed parameters (we are using the simplest form of the 
simulated Y500 — M500 relation here). 

Another interesting result from the PXMM sample is 
that the results are consistent with Ysoo = Yx.soo on aver- 
age (again like our PC and FO models), however the scat- 
ter in the observed relation is significantly larger than ours 
(observationally, o"i O g 10 Y = 0.1, around a factor of 5 larger 
than for our PC and FO simulations). As a result, the scat- 
ter in the other observed PXMM scaling relations are also 
larger than ours; e.g. the scatter is 2-3 times larger for the 
Y500 — M500 relation. Thus if our PC and FO simulations, 
calibrated to X-ray data, provide faithful estimates of the 
mean SZ/X-ray scaling relations, observational estimates of 
the quantities must somehow increase the scatter without 
introducing significant bias. One potential source of scatter 
is due to the projection of large-scale structure along the 
line-of-sight; we investigate this below. 



4.6 Projection effects 

As detailed in Section 12.61 we have constructed 50 5° x 5° 
mock realisations of the SZ sky (Compton y maps) from our 
GO and PC simulations. (Unfortunately, it is not currently 
possible to do this for the FO model as it was not run on 
the full Millennium volume.) We use these maps to estimate 
the (cylindrical) Y500 for the clusters that are present, as 
follows. 

Firstly, we cross-match our 50 maps with cluster cata- 
logues at all available redshifts (catalogues are constructed 
for all snapshots used to make the maps, providing there are 
objects above our mass limit of M500 = 10 14 Ii~ 1 Mq). This 
is done by performing the same operations (translation, ro- 
tation, reflection) on the cluster centre co-ordinates as was 
done with each of the snapshots, then finding the pixel in 
the map that corresponds with the cluster centre, for those 
objects within the map region. We then identify which pixels 
fall within the projected radius, -R500 = ?"500, and compute 
the SZ Y500 value which we define as 



Y. 



sky 



D% 8Q. ^ Vi,i 



(18) 



where the sum is performed over all relevant pixels (with 
indices, i,j) and SSI is the solid angle of each pixel (we use 
1200 x 1200 pixels so assume Sfl = 0.25 x 0.25 arcmin 2 ). 
Finally, we throw away clusters that have a more massive 
neighbour whose centre lies within its own radius, -R500, 
as this interloper would dominate the estimated SZ flux. 
Our final catalogue is restricted to clusters with M500 > 
10 14 /i _1 M0 and z < 1; for comparative purposes we split 
this into a low-redshift (z < 0.5) and high-redshift (z > 0.5) 
sub-samples. The number of clusters in each of these sub- 
samples for the GO and PC models are listed in Table|4] The 
larger numbers in the high-redshift sample are expected due 
to the larger volume there (for fixed solid angle). Note that 
the same cluster could appear more than once (in a different 
realisation or redshift). 

In order to extract the cluster signal from the rest of the 
large-scale structure along the line-of-sight, we also compute 
cylindrical Y500 values due to the cluster region itself. To do 
this, we apply equation (|18[1 to our cluster maps, detailed 
in Section 12.51 As a reminder, the length of the cylinder, 
centred on the cluster, is z = 12r5oo; this approximately 
corresponds to three virial radii from the centre in each di- 
rection along the line-of-sight. We refer to this Y value as 
Y 5 cl r; clearly Y 5 s k y > Y^ s by definition. 

The squares in Fig. [15] represent the YsocT — Afsoo rela- 
tion for our GO (top panels) and PC (bottom panels) mod- 
els at high (left panels) and low (right panels) redshifts. We 
re-scale cluster Ysoo s values by E(z)~ 2 ^ 3 to account for evo- 
lution across the redshift range in each panel. Best-fit pa- 
rameters (A, B, criogY) are given in Table [4] a pivot mass of 
3 x 10 14 /i _1 M Q was adopted for all the fits. The GO model 
relations show similar trends to those seen in the spherical 
Y500 — M50Q relation; the slope is close to self-similar and 
the scatter is small. The PC relations again have slopes that 
are steeper than the self-similar value but also have slightly 
larger scatter (o"i ogl0 y — 0.07 — 0.08), reflecting in part the 
effect of additional evolution with redshift. 

The stars in Fig. [15] are for when Y£ q values are used 
and thus contain the additional signal from beyond the clus- 
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Table 4. Best-fit parameters for simulated SZ scaling relations using projected (cylinder) values from cluster and sky maps (see text for 
further details). Column 1 gives the scaling relation being considered; column 2 the redshift range; column 3 the simulation model; column 
4 the number of clusters used in the fit; columns 5 & 6 list the best-fit values for the normalisation and slope parameters respectively; 
and column 7 lists the estimated scatter in o"i oglQ y . Quoted uncertainties correspond to either the 16th or 84th per centile (whichever 
is largest), estimated using the bootstrap re-sampling technique. 



Flux 




Redshift 


Model 


N i 

1 v clus 


A 


13 


°"logio Y 


v^clus 
J 500 


M 500 


< z < 0.5 


GO 


1346 


-4.677 ± 0.003 


1.650 ±0.007 


0.045 ± 0.001 






PC 


1074 


-4.671 ± 0.004 


1.72 ±0.01 


0.068 ± 0.002 






0.5 < z < 1 


GO 


2952 


-4.702 ± 0.002 


1.613 ±0.007 


0.050 ±0.001 








PC 


2199 


-4.699 ± 0.003 


1.74 ±0.01 


0.077 ±0.001 


v sk y 

-^500 


M 500 


< z < 0.5 


GO 


1346 


-4.622 ± 0.003 


1.507 ±0.009 


0.059 ±0.001 






PC 


1074 


-4.455 ± 0.003 


1.23 ±0.01 


0.059 ±0.001 






0.5 < z < 1 


GO 


2952 


-4.677 ± 0.003 


1.524 ±0.007 


0.062 ±0.001 








PC 


2199 


-4.556 ± 0.003 


1.32 ±0.02 


0.071 ±0.001 
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r 500 
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< 2 < 0.5 


GO 


1346 


-4.686 ± 0.004 


1.71 ±0.01 


0.074 ± 0.003 






PC 


1074 


-4.683 ± 0.005 


1.80 ±0.02 


0.106 ±0.003 






0.5 < z < 1 


GO 


2952 


-4.719 ± 0.003 


1.666 ±0.009 


0.071 ±0.001 








PC 


2199 


-4.721 ± 0.004 


1.80 ±0.01 


0.104 ±0.002 




1 10 1 10 1 10 1 10 

Msoo/1 O'V'M,, M 500 /1 , 4 h _, M o M 500 /1 , 4 h _, M o M 500 /1 , 4 h _, M o 



Figure 15. Projected I500— relations for clusters in the GO 

(top panels) and PC (bottom panels) sky maps with < z < 0.5 
(left panels) and 0.5 < 2 < 1 (right panels). The stars correspond 
to values, i.e. calculated from the full sky map. The squares 

correspond to Kfgg 8 values, i.e. from the cluster region. In both 
cases the true M500 values were used. The dashed line is a best- 
fit to the I^qq 8 — M500 relation and the solid line to the Y^ q — 
M500 relation. The Y^qq values are higher on average than Vggg 8 , 
especially in the PC simulation at low mass and low redshift, 
where the difference is a factor of 2-3. 

ter. The difference between the two relations in each panel 
(as can be seen from the best-fit lines) is most prominent 
for the PC model, where the slope has decreased from ~ 1.7 
to ~ 1.2, due to Y^qq being significantly larger than Y^o 8 
in the lower mass objects. As was discussed in Section [2] 
the pre-heating was applied everywhere at z — 4 and thus 



Figure 16. As in Fig. I15l but the Y^ Q q values have had the mean 
background signal subtracted. The two best-fit relations are now 
very similar in all panels. 



substantially increased the thermal energy of the gas, as in- 
dicated by the three-fold increase in the mean y signal. Such 
widespread heating is likely to be unrealistic as it would 
require a huge amount of energy and would boil off the 
small amount of neutral hydrogen and helium in t he IGM 
l|Theuns, Mo fc Schavell200ll ; iBorgani fc Vielll2009l ). so the 
PC result represents a worse-case scenario for the effects of 
projection on the Y signal. 

Observations of the SZ effect made with the Planck 
satellite are unable to measure the mean y signal as at 
each frequency, spatial temperature fluctuations are mea- 
sured with respect to the all-sky mean. It is therefore more 
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realistic to compare the background-subtracted values of 
to the cluster values. To do this we compute the projected 
angular area for each cluster and compute the expected con- 
tribution to Y500 from the mean y 



y^skysub _ ysky _ ^ ^2 



(19) 



where fisoo is the solid angle subtended by the cluster out to 
a projected radius, ifeoo- The results of this procedure are 
shown in Fig. 1161 with best-fit parameters for the igoo"" 1 — 
A/500 relations given in Table 2] 

Interestingly, the two best-fit relations are now almost 
identical for each run and within each redshift range. A sim- 
ple background subtraction therefore removes any bias in the 
mean relation generated from the additional hot gas along 
the line-of-sight. The scatter is considerably larger in the 
Y 5 s k ysub — M500 relation, in part due to the fact that the ad- 
ditional signal is not constant everywhere. The PC relations 
again contain the largest scatter, comparable to the observed 
scatter in the PXMM data (<xi ogl0 y ~ 0.1). Although the re- 
sult is model dependent, it is clear that part (if not all) of 
the observed scatter can be attributed to projection effects. 



5 HYDROSTATIC BIAS 

In the previous section, we saw that our PC and FO mod- 
els produced SZ/X-ray scaling relations that were in good 
agreement with the PXMM observational data. A significant 
uncertainty in the observational determination of scaling re- 
lations is the (direct or indirect) assumption of hydrostatic 
equilibrium (HSE), required for deriving the cluster mass 
(•W500) and radius (rsoo). It is therefore interesting to look 
at the accuracy of this assumption in our simulations as the 
good agreement between our results and the observations 
can only be preserved if hydrostatic bias is small (in the 
absence of additional systematic effects). 

For a cluster in HSE, the pressure gradient in the ICM 
is sufficient to balance gravity; the total mass of the cluster 
can then be calculated as 



M 



\< r) = - 



kTr 
GfrniH 



dlnp dlnT 

+ 



d In r d In r 



(20) 



where /1 = 0.59 is the mean molecular weight for an 
ionised plasma (assuming zero metallicity) . We use the 
spectroscopic-like temperature to evaluate the local temper- 
ature, T(r) and its gradient, dlnT/dlnr, at radius, r. 

Estimation of the cluster mass based on hydrostatic 
equilibrium can be biased for three reasons. Firstly, the es- 
timated mass within a fixed radius can be different from the 
true mass because the intracluster gas is not perfectly hydro- 
static. Previous simulations have shown that mass estimates 
can be too low by up to 20 per c ent, due to incomplete ther- 



malisation 
19961; 



of the gas (e.g. 



Kay et ah! 
20071; 



IRasia. Tormen fc Moscardin: 



I2004 



Evrard. Metzler fc Navarrd 
200? 



IRasia et al 



2006 



Nagai, Vikhlinin & Kravtso- 



Nagai Kravtsov fc Vikhlininl 120071; iPiffaretti fc Valdarninil 



iKav et all 
20071; 



20081 ; Umeglio et all 120091 ; iLau. Kravtsov fc Nagail 12009) . 
A second effect is that the X-ray temperature of the gas 
may be lower than the mean (mass-weighted) temperature. 
Such an effect depends on the thermal structure of the 
gas (in particular, the low entropy tail associated with 



substructure) and can be particularly severe when radiative 
cooling effects are strong, [f] Finally, the cluster's size itself 
is usually defined as a scale radius (e.g. rsoo) which is 
mass-dependent so also depends on the assumption of 
hydrostatic equilibrium. 

To study how these combined effects impact upon our 
scaling relations, we estimate the hydrostatic mass of each 
cluster as follows. Firstly, we compute the hot gas (T > 
10 s K) density and temperature profiles. In lower mass clus- 
ters the profiles can get rather noisy due to limited particle 
numbers which can affect the estimation of the pressure gra- 
dient. To avoid this, we fit a cubic polynomial function to 
each profile (in log space) to generate a smoothed represen- 
tation. (This also has the advantage that the gradient can 
be derived analytically.) We then use these model profiles 
to estimate the mass, M HSE , using equation (|20[) . then vary 
the radius, r, until the following equation is satisfied 



Mr, 



4tt 



(rg?, E ) 500^(2), 



(21) 



where M|J E and r^§^ are our estimated mass and radius 
respectively. Once the radius is known we can use this to es- 
timate the SZ flux which we will denote Y^ E . A gain, this is 
the flux from within a sphere centred on the cluster; all that 
has changed is the assumed value of rsoo - In what follows, we 
only consider the sub-set of clusters in the estimated mass 
range, lO 14 ft _1 M < Mjo E < 1O 15 /i _1 M . The numbers 
of clusters are listed for each model and redshift in Table [5] 



5.1 Effect of hydrostatic assumption on cluster 
mass 

We quantify the effect of the hydrostatic assumption on clus- 
ter mass by considering the distribution of the estimated-to- 
true mass ratio, Rm = log 10 (M|J, E /Af50o), for our models 
at z = l and z = 0. (Note that _Rm directly measures the 
resulting shift along the logarithmic mass axis.) The results 
are shown in Fig. 1171 

At z = 0, the GO results show a significant spread in 
mass ratios as well as a large negative bias; the median value 
is _Rm = —0.14. In the PC and FO models, the spread and 
bias is smaller, with the median increasing to around —0.05. 
A similar situation is evident at z = 1. The disturbed sub- 
sample, where HSE should definitely not be a good approx- 
imation, shows a small offset in the median Rm from the 
overall sample; in the PC and FO cases the offset is positive 
whereas in the GO case it is negative. 

It is perhaps not surprising that the discrepancy be- 
tween estimated and true mass from the GO simulation is 
significantly higher than for the PC and FO models. As is 
evident from the Y500 — Yx,soo relation (Fig. ll3|l . the former 
model predicts a more clumpy intracluster medium due to 
the persistence of low entropy gas that is unable to cool. This 
gas by its very nature has not completely thermalised to the 
global cluster temperature and has significant residual bulk 
kinetic energy. In the latter two runs, the non-gravitational 



6 We note that lNaeai. Vikhlinin fc Kravtsovl j2007h found the X- 
ray temperature to be higher than the mass- weighted temeprature 
in a mock Chandra analysis of their simulated clusters, but they 
exclude any resolved cold clumps from their calculation. 
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Table 5. Best-fit parameters for simulated SZ scaling relations, with estimated properties from the assumption of hydrostatic equilibrium. 
Column 1 gives the scaling relation being considered; column 2 the redshift; column 3 the simulation model; column 4 the number of 
clusters used in the fit; columns 5 &; 6 list the best-fit values for the normalisation and slope parameters respectively; and column 7 lists 
the estimated scatter in &i o g 10 Y or CT log 10 M (whichever is appropriate). Quoted uncertainties correspond to either the 16th or 84th per 
centile (whichever is largest), estimated using the bootstrap re-sampling technique. 



Relation 




Redshift 


Model 


N i 

1 ' cms 


.4 


B 


°"log 10 Y/M 




- M HSE 

1V1 500 


z = 


GO 
PC 
FO 
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179 


-4.54 ±0.02 
-4.75 ± 0.009 
-4.69 ±0.01 


1.60 ±0.07 
1.69 ±0.03 
1.50 ±0.05 
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0.44 ±0.01 
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0.58 ±0.02 
0.545 ±0.009 
0.60 ±0.02 


0.118 ±0.004 
0.064 ± 0.003 
0.080 ± 0.007 






2 = 1 


GO 
PC 
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0.37 ± 0.06 
0.69 ± 0.03 
0.56 ± 0.02 
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0.57 ±0.04 
0.59 ±0.04 
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0.066 ± 0.008 
0.09 ±0.01 
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-4.356 ± 0.009 
-4.694 ± 0.001 
-4.668 ± 0.002 


1.84 ±0.02 
1.840 ±0.004 
1.659 ±0.008 


0.055 ±0.002 
0.0150 ± 0.0005 
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2 = 1 


GO 
PC 
FO 


31 

102 
75 


-4.29 ±0.03 
-5.040 ± 0.002 
-4.799 ±0.004 


2.09 ± 0.09 
1.757 ±0.005 
1.67 ±0.02 


0.052 ±0.007 
0.0092 ± 0.0008 
0.022 ±0.006 



heating generates a smoother distribution that is evidently 
closer to hydrostatic equilibrium. 



5.2 Effect of hydrostatic assumption on Y500 and 

Yx,500 

The use of hydrostatic mass estimates also affects the SZ 
flux through the use of r|oo E to define the cluster radius; a 
smaller radius will result in a lower value for Y. We define a 
similar quantity to the mass ratio, Ry = log 10 (Y 5 oo E /Y5oo), 
and present the distribution of values in Fig. 1181 Again, we 
present the ratio in this way as it directly gives the shift in 
log 10 Y values due to the hydrostatic estimate. 

As was the case with the total mass estimates, there 
is a larger bias (and scatter) in the Y500 values for the GO 
run but the overall effect is smaller as it is entirely due to 
the (small) shift in rsoo- The median Ry is —0.03 for GO 
at z = 0, increasing to only —0.02 for the PC and FO runs. 
Since r^ o E < ^500 on average, the integrated flux is also 
smaller. Again, the results are not significantly different at 
high redshift or when only the disturbed clusters are se- 
lected. We have also checked the equivalent result for the 
Yx,500 values and they are very similar to the Y500 results. 



5.3 Estimated Y500 — Afsoo relation directly from 
HSE 

We now put together these results to study how the Y500 — 
M500 relation is affected by the hydrostatic assumption. 
These results are shown in Fig. [19] for the GO, PC and 
FO runs at z = 1 and z = 0. The best-fit Y^ - M|^ E 
relation is shown as a solid line and we also plot the best-fit 
(true) Y500 — M500 relation as the dashed line in each panel. 



Values for the parameters describing the best-fit relations 
(normalisation, A; slope, B; and scatter, o"i ogl0 Y") are given 
in Table [5] 

The offset in M500 values (Rm = —0.14 at z = 0) in the 
GO model is clearly visible in the top-right panel of Fig. 1191 
where the best-fit relation is offset to larger Y500 values for a 
given value of Afjjoo • The large spread in the Rm distribu- 
tion is also evident as the scatter has increased significantly 
(cri ogl0 Y ^ 0.19, c.f. Fig.|S]where <7 logl0 Y ^ 0.04). The offset 
is insensitive to mass in this model, resulting in a relation 
that has similar slope (1.6) to the true Y500 — M500 relation. 
The offset in normalisation has also led to a significant drop 
in the number of clusters in the sample at each redshift; as 
a result there are only 98 clusters at z = 1, making a re- 
liable estimate of the relation difficult (but the trends are 
nevertheless consistent with those seen at 2 = 0). 

The best-fit Y^ E - relation from the PC run 

at z = is remarkably similar to the underlying rela- 
tion, although the scatter has also increased considerably 
to criog 10 y — 0.11. Results at z = 1 prefer a flatter slope but 
this is somewhat affected by a few higher mass clusters (the 
slope is 1.5 ±0.1). The estimated relation for the FO model 
is also similar to the true relation, with a preference for a 
slightly flatter slope and larger scatter (ai ogl0 y ~ 0.13 at 
2 = 0). The disturbed cluster sub-sample is most strongly 
biased in the PC results at z = 0, where the clusters have 
larger HSE masses for their flux, relative to the regular sys- 
tems. 

5.4 Estimated V500 — M500 relation using Ix.soo 

When mass estimates are required for larger samples of 
clusters, the direct hydrostatic method discussed above 
can be prohibitively expensive as it requires the den- 
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Figure 17. Distribution of estimated-to-true mass ratios, i?M; 
within the estimated rsoo for clusters at z = 1 (left panels) and 
2 = (right panels). The top panels show results from the GO 
model, middle panels from PC and bottom panels from FO. The 
green histogram is for the whole cluster sample while the blue 
histogram is for the disturbed sub-sample. Vertical dashed lines 
indicate the median mass ratio for each case, with values given 
in the legend. The median Hm is significantly smaller in the runs 
with non-gravitational heating. 



sity and temperature profiles to be known out to r 5 oo 
and beyond. An alternative, indirect method is to use 
a mass proxy, where mass is estimated from a mass- 
observable scaling relation that is pre-calibrated using 
fewer clusters. Historically, Tx was the observable of choice 
but recent studies have focussed on the use of Yx due 
to i ts low scatter dKravtsov. Vikhlinin fc Nagail 120061; see 



also I Arnaud. Pointecouteau fc PrattJ 120071; iMaughanl 120071: 



lArnaud et all 120101: ISun et all 120111). Indeed, the Planck 



Collaboration ( Planck Collaborationll2pild) made us e of the 
Msqo — Yx,soo relation, calibrated by Arnaud et alj (|2010l ) 
from the REXCESS sample of 33 clusters, to estimate rsoo 
and M500 for their larger (PXMM) sample of 62 clusters. 

The procedure for estimating M500 works as follows. 
Assuming that all clusters lie on an M500 — Yx, 500 relation 
and that they evolve self-similarly with redshift, then rsoo 
may be found using 

rsooCKx) = a E{ Z y i/5 m(F x ) 1/3 ft^Mpc, (22) 



where a is a known constant and m(Yx) is the best-fit value 
of E 2 ^ 5 Mc,oo(Yx,50q) from the scaling relation for a given 
Yx,500- The mass, Msoo(5 / x), can then be estimated using 
equation (|21[) . In practice, this equation must be solved it- 



Figure 18. As in Fig. 1171 but for the ratio of estimated-to-true 
^500 values, Ry- The difference in log(Ysoo) is very small in the 
runs with non-gravitational heating. 



eratively: a value for rsoo is first guessed then Yk,soo is cal- 
culated within this radius (from the integrated gas mass and 
average X-ray temperature), allowing a new value for rsoo 
to be computed from equation (1221) . This is repeated until 
convergence is achieved. Since clusters do not all lie on this 
relation (even though this particular relation is chosen for 
its low scatter) the derived rsoo may be inaccurate for an 
individual cluster, but the overall relation should be unbi- 
ased. 

We have applied this procedure to our simulated clus- 
ters and will refer to the resulting Y500 — M500 relation 
as the Yg^Q — relation. We first show our derived 

M^o E - ix^500 relations, required for equation (f22|) . for the 
three models at z = 1 and z = in Fig. [20] Since Yx (and 
not mass) is on the :r-axis, we restrict our fits to clusters 
with 3 x 10 13 /i~ 1 M Q keV < Y^ < 10 15 /i" 1 M keV, 
as this approximately matches our adopted mass range 
(10 14 - 10 15 h~ 1 M Q ) for the PC and FO models. 

Qualitatively, the same conclusions can be drawn as 
for the Y£o E - M^ E relation (Fig. OH): the GO model 
shows a large scatter and the relation is offset due to the 
mass estimates being systematically low. As expected, how- 
ever, the PC and FO results agree very well with the best- 
fit underlying relation with still relatively small scatter, 
ciog 10 m = 0.06 — 0.08. We also compare our results to the 
observ ed best-fit relation at low redshift from I Arnaud et all 
(2010), shown as the dot-dashed line. The agreement be- 
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Figure 19. The T 5 qq E — *^loo E relations (i.e. using estimated 
values for each cluster assuming the gas is hydrostatic) for the 
GO, PC and FO models at z = 1 and z = 0. The squares corre- 
spond to regular clusters while the triangles are disturbed clus- 
ters. The solid line is the best-fit relation while the dashed line 
shows the best-fit true Y500 — M500 relation. Fits are performed 
for all clusters with 10 14 A^Mq < < 10 15 A^Mo and 

best-fit parameter values are given in Table [5] It is clear that the 
hydrostatic assumption is more robust for the PC and FO runs 
than for the GO run. 

tween the PC and FO models is very good, with the latter 
prefering a slightly flatter slope than the observed relation. 

The derived Y^j — relations are shown in Fig. 1211 
We have overlaid the best-fit relation to these data (solid 
line); the best-fit true Y500 — A^soo relation (dashed line); 
the best-fit y 5 oo E - M^ % E relation (triple-dot-dashed line); 
and the best-fit relation to the PXMM observational data 
(dot-dashed line). 

The final result is very striking for our PC and FO mod- 
els. At z = 0, there is only a small amount of bias (around 
20 per cent or so, or A log 10 Y500 — 0.08) with respect to the 
underlying relation (and also, the PXMM relation) . Such an 
offset, comparable to the estimated intrinsic scatter in the 
observed relation, is small enough that it may just reflect our 
method not exactly matching that used by the Planck Col- 
laboration. Even more striking is the reduction in scatter; 
°"iog 10 y — 0.02 for the PC and FO models at z — 0, about 
half the size of the scatter in the true Y500 — M500 rela- 
tion. The reason for this reduction is obvious: from Fig. [13] 
we saw that Y500 and Ix.soo were strongly correlated, es- 
pecially in the PC and FO models. Thus, a cluster that 
has a larger-than-average Yx.soo for its mass will also have 



Figure 20. As in Fig. [19] but for the M|^ E - Y^| E relation. 

Clusters are now selected with 3 x 10 13 < Yxsoo < 1Ql5 wnicn > 
for the PC and FO models, matches well to our normal mass 
range. The solid line is the best-fit to the relation and the dashed 
line the true relation. The dot-dashed line is the b est-fit relation 
from R EXCESS (XMM-Newton) data as found bv lArnaud et al.l 

teoich . 

a larger-than-average Igoo- This was also true for the dis- 
turbed sub-sample hence the reason why these clusters are 
unbiased with respect to the overall sample. 



6 SUMMARY AND CONCLUSIONS 

Large surveys are now being performed at millimetre wave- 
lengths exploiting the Sunyaev-Zel'dovich (SZ) effect to de- 
tect large samples of galaxy clusters out to high redshift. 
Such samples will then be used to produce competitive con- 
straints on cosmological parameters, as well as to study the 
variation in physical properties of the intracluster gas (es- 
pecially the gas pressure) with mass and redshift. The cos- 
mological application relies on the statistical estimation of 
clus ter mass through the S ZY — M relation. In recent w ork 
(e.g. lAndersson et al.|[201ll ; |Planck CoUaborationll2011ch the 
first SZ-selected samples of clusters have already been used 
to estimate the cluster Y — M relation and full results of 
cosmological analyses are expected over the next few years. 

In this paper, we have analysed some of the largest N- 
body/hydrodynamic simulations of structure formation (the 
Millennium Gas Simulations) to study the dependence of SZ 
cluster properties on gas physics, at both low [z = 0) and 
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Figure 21. Y500 — M500 relations for the three models at z = 
1 and z = when rsoo/jWsoo is estimated from the best-fit 
-^500 E — ^x^soo relation with no scatter. We also show the best-fit 
true Y500-M500 relation (dashed line), the best-fit Y^ E -Mf^ 
relation (triple-dot-d ashed line) and the best-fi t 2 < 0.5 relation 
from PXMM data l|PIanck Collaborationll201ld : dot-dashed line). 
The PC and FO models agree well with the observational data al- 
though there is a slight offset. The scatter is much lower, however, 
due to the strong correlation between Y500 and yx,500- 



high (z — 1) redshift. The large volume used in these sim- 
ulations produces significant (hundreds to thousands) sam- 
ples of clusters over the interesting range of cluster masses 
(10 14 — 10 15 /i _1 M0). We considered three cluster gas physics 
models: a non-radiative (gravitational heating only) simu- 
lation that ought to produce an approximately self-similar 
cluster population; and two simulations that incorporate ad- 
ditional non-gravitational heating (a model that uniformly 
pre-heats the gas at high redshift and a model that includes 
feedback from stars and active galactic nuclei in galaxies). 
The feedback model is our most realistic, in that it has al- 
ready been shown to reproduce many of the scaling prop- 
erties of X-ray clu sters, especially those with non-cool cores 
jShort et alj|201Ch . 

We started by investigating the hot gas pressure pro- 
files of our simulated clusters and how they comp are to the 
pressure profile advocated bv lArnaud et al.l |2010l ). We then 
compared our derived SZ scaling relations (between Y500 
and total mass, hot gas mass, X-ray temperature and the 
X-ray analogue to the SZ Y parameter, Yx) with the recent 
observational results, in particular those obtained from a 
combined SZ+X-ray analysis performed by the Planck Col- 



laboration. We also tested two of the key assumptions used 
in the observed analysis, namely that the mean I500 — A/500 
relation is unaffected by the assumption that the gas is hy- 
drostatic and by the presence of any other hot gas along the 
line-of-sight. Our main conclusions can be summarised as 
follows: 

• In accord with previous studies, our simulation with 
non-radiative hydrodynamics produces a (spherical) F500 — 
M500 relation that has a self-similar slope (5/3) and also 
evolves with redshift according to the self-similar expec- 
tation, E(z) 2 ^ 3 . Simulations with non-gravitational heating 
(both pre-heating and feedback cases) create slightly steeper 
Y500 — M500 relations (with a slope of 1.7—1.8, when clusters 
across the mass range, 10 14 h~ 1 M Q < M 500 < 10 15 h~ 1 M Q , 
are considered) but the evolution with redshift is still close 
to self-similar. 

• The simulations were compared with the Pla nck+XMM 
results at z < 0.5 (|Planck Collaboration! l201ld ) and very 
good agreement was found for a number of scaling relations 
(isoo versus M500, Af gaa ,5oo and kTx) for the pre-heating and 
feedback models. The scatter in the Y500 — M500 relation is 
smaller than observed, however, with o"i ogl0 y — 0.04. 

• Intracluster gas in the non-radiative simulation contains 
a significant unthermalised component, due to the presence 
of low-entropy, clumpy gas. This causes an offset in the 
Y500 — ix,soo relation, which tests the difference between 
the mass- weighted and X-ray temperatures. As a result, hy- 
drostatic mass estimates are biased low by 20-30 per cent. 
The pre-heating and feedback simulations on the other hand 
predict smoother gas distributions, with Y500 — Yx,soo and 
much smaller hydrostatic bias (estimated masses are only 
~ 10 per cent lower). 

• The estimated M500 — Yx,soo relations (assuming the 
gas is hydrostatic) are in goo d agreement w i th th e recent 
observational determination bv lArnaud et al] (|2010T) . When 
^x,500 is used as a mass proxy to predict the SZ Y500 — M500 
relation, only a small (~ 20 per cent) offset in normalisation 
from the true relation (and thus the observed relation from 
Planck+XMM data) is found. The scatter in the recovered 
relation is very small (ciog 10 Y — 0.02) due to the strong 
correlation between Y500 and Yx, 500- Clusters that are un- 
dergoing major mergers are not significantly offset from the 
mean relation. 

• Hot gas pressure profiles are well described 
by generalised NFW p r ofiles, as suggested by 
iNagai, Kravtsov fc Vikhlininl l|2007h and show that the 
majority of the contribution to the SZ Y parameter (where 
r 3 P(r) is maximal) comes from radii close to rsoo. Splitting 
the cluster samples into low and high-mass sub-samples, 
we find little difference between the two in the run with 
non-radiative hydrodynamics, as expected. The runs with 
non-gravitational heating predict that low-mass clusters 
have lower core pressures and higher pressures in the cluster 
outskirts, when scaled according to the self-similar expec- 
tation. This non-self-similar-behaviour can be attributed 
to the heating that is more effective in low-mass clusters 
and acts to push the gas out to large radii. There is also 
significant cluster-cluster scatter, especially in the core 
region and in the outskirts, where individiual pressure 
profiles can be 50 per cent higher than the median profile. 

• We also compared our median pressure profiles with 
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the Arnaud et al. profile and found good agreement (within 
10 per cent) for our high- mass clusters at r > 0.5rsoo, in 
the pre-heating and feedback models. Low-mass clusters are 
especially discrepant in the core regions, likely due to the 
absence of radiative cooling in our models. Using the X-ray 
temperature (rather than hot gas mass-weighted tempera- 
ture) in the pressure calculation, as well as using hydrostatic 
estimates of rsoo and M500, only makes a significant (> 10 
per cent) difference to the non-radiative simulation for the 
reasons already mentioned. 

• Finally, we considered the effects of projection due to 
large-scale structure along the line-of-sight, by analysing 50 
5° x 5° maps of the thermal SZ effect. By measuring the 
cylindrical SZ flux associated with each cluster and com- 
paring to the flux from the cluster region alone, we were 
able to discern the contribution from additional structures, 
in the non-radiative and pre-heating simulations. The pre- 
heating model showed the largest bias, where low-mass clus- 
ters (M500 — 10 14 /i _1 Mq) had cylindrical Y500 values that 
were around 2-3 times higher than the value from the cluster 
region. This is due to the large amount of thermal energy 
injected into the gas at high redshift, as evidenced by the 
three-fold increase in the mean-y parameter. Subtracting the 
contribution from an assumed mean background we find the 
recovered K500 — M500 relation to be unbiased with respect 
to the cluster relation, with some additional scatter that is 
mo del- dependent. 

In summary, we can conclude that when our more re- 
alistic models for the intracluster gas are employed (namely 
those that raise the entropy of the gas to match global X-ray 
scaling relations), the SZ Y — Al relation is in good agree- 
ment with the observations (Fig. [Sj and is largely unaffected 
by two of the main sources of systematic uncertainty: hy- 
drostatic bias (Figs. [T9"l and |2"T| and projection effects from 
large-scale structure (Fig. [T6)l . 

While our analysis has been one of the most comprehen- 
sive to date and used some of the largest and most sophis- 
ticated simulation models, there are some significant short- 
comings that still need to be addressed. Firstly, the effects 
of radiative cooling were not included in our most realistic 
(feedback) model, so the model cannot yet match the full 
X-ray clust er population (na mely the brightest objects with 
cool cores. Ishort et aT1l2010h . As we argued, this omission 
is likely not a significant problem for the Y — M relation 
but will affect the hot gas pressure profile so it should be 
addressed in future work. Secondly, we were unable to test 
projection effects for the feedback model as we only have 
a sample of clusters rather than the full cosmological vol- 
ume. Finally, the cosmological model adopted for the sim- 
ulations (identical to that used in the original Millennium 
Simulation) is no longer favoured; in particular the value 
of as is higher than the current best estimate (as = 0.9 in 
the simulations, c.f. as ~ 0.8 from the WMAP 7-year data; 
iKomatsu et al.l 201lh . Using the presently- favoured cosmo- 
logical model is likely to reduce the scale of projection ef- 
fects, however, as in it structure formation will be less ad- 
vanced. 

We are currently preparing a new generation of Mil- 
lennium Gas simulations that will rectify all of these prob- 
lems, starting with a new version of our existing feedback 
model that will deal with the second and third issues. This 



new simulation, which is also being run at higher resolution 
and with an upda ted semi-analytic galaxy formation model 
|Guo et al.l l2010f ). will additionally allow cosmologically- 
dependent statistical predictions for the SZ signal to be per- 
formed, namely the SZ power spectrum. 
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Figure Al. As in Fig. 1101 but for the Y500 — A^gas,500 relation. 



APPENDIX A: EVOLUTION OF SCALING 
RELATIONS 

The following figures (Figs. [Al] |A2] and |A3j illustrate the 
evolution of the slope, normalisation and scatter with red- 
shift for the y 5 oo - M gaSi5 oo, Y 500 - T A and Y 50 o - Vx.soo re- 
lations respectively. Details of what is plotted in each panel 
are identical to Fig. [TU] and are discussed in Section 14.21 

This paper has been typeset from a TpX/ DTgX file prepared 
by the author. 
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Figure A3. As in Fig. UOI but for the Ygoo — 500 relation. 
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